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CHAPTER  1 
INTRODUCTION 

1.1  Overview 

This  thesis  is  an  investigation  into  the  use  of  different  types  of  Global 
Positioning  System  (GPS)  data  together  with  different  types  of  bias  sources.  The 
intent  was  to  examine  their  relative  effect  on  the  determination  of  the  position  of 
the  receiver  and  the  accuracy  of  the  retrievable  information.  Selected  inputs  were 
processed  through  a  simulation  software  package,  "Multimode  Simulation  for 
Optimal  Filter  Estimation  (MSOFE)”  (Carlson  and  Musick,  1990).  GPS  users 
and  researchers  continually  attempt  to  analyze  the  effects  of  data  and 
measurement  biases  to  increase  the  quality  of  information  gained  from  the 
system.  In  this  thesis,  pseudorange  and  phase  range  data,  together  with  various 
estimated  and  ignored  measurement  biases,  were  simulated  and  then  processed 
through  the  optimal  estimation  software  package.  The  output  of  the  software- 
optimal  estimations  of  the  desired  quantities-provided  insight  into  the 
consequence  of  the  virtually  endless  combinations  of  data. 


» 


1.2  Basic  Statement  of  Problem 


Figure  1  is  a  graphic  representation  of  the  underlying  principle  behind  this 
thesis.  What  happens  to  our  ability  to  recover  positions  when  the  most-simple, 
ideal  situation  is  augmented  with  other  forms  of  data  and  perturbed  by  biases? 
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FIGURE  1 

PRINCIPLE  BEHIND  THESIS 

The  preliminary  stage  utilized  only  the  most  basic  GPS  measurement  type, 
the  pseudorange,  which  is  known  to  about  one  meter.  The  second  stage,  the 
incorporation  of  the  phase  range,  known  to  about  a  millimeter,  attempted  to 


increase  the  accuracy  of  the  results.  Minus  any  biases  or  errors,  these 
combinations  comprised  the  two  "ideal"  cases.  To  these  ideal  cases, 
perturbations  in  the  form  of  biases1  will  be  added  piecemeal.  The  biases  will  fall 
into  three  different  categories: 

a.  station-dependent  biases— i.e.,  station  coordinates  in  x,  y,  and 
z;  and  receiver  clock  biases2; 

b.  observation-dependent  biases— i.e.,  tropospheric  delay,  iono¬ 
spheric  delay,  and  phase  integer  ambiguity  biases;  and 

c.  satellite-dependent  biases— i.e.,  satellite  clock  and  radial  orbit 
biases. 

The  study  of  adding  additional  data  and  biases  will  be  accomplished  with 
the  MSOFE  software  package.  This  software  provides  for  Kalman  filtering-a 
sequential,  recursive,  mathematical  data  adjustment  procedure  designed  to 
provide  optimal  estimation  of  errors  within  a  linear  system.  Gelb  (1992) 
describes  an  optimal  estimator  as  that  which  is  generated  by 

a  computational  algorithm  that  processes  measurements  to  deduce  a 
minimum  error  estimate  of  the  state  of  a  system  by  utilizing:  knowledge 
of  system  and  measurement  dynamics,  assumed  statistics  of  system  noises 
and  measurement  errors,  and  initial  condition  information. 


]In  this  thesis,  the  term  "bias'  will  be  used  for  those  random  effects  on  the  measurement  which 
must  be  modeled  (systematically)  and  corrected  for  in  order  to  obtain  an  accurate  GPS  range 
measurement  (Wells,  1987). 

2A!though  these  quantities  are  typically  estimated  in  a  GPS  adjustment,  they  are  considered  biases 
here  because  they  influence  the  GPS  measurement  and  must  be  accounted  for  as  such. 


Gelb  further  states  that  this  form  of  estimation  has  some  distinct  advantages. 
First,  the  data  processor  (here,  the  computer  program)  is  able  to  minimize  the 
estimated  error  in  a  well-defined,  statistical  sense.  Second,  the  processor  is  able 
to  utilize  all  measurement  data  along  with  any  a  priori  knowledge  about  the 
system  (Gelb,  1992).  This  estimator  is  said  to  be  used  in  a  "filtering"  capacity 
when  the  time  at  which  the  information  is  desired  is  the  last  available 
measurement  epoch.  That  is,  no  future  data  are  available  for  use  in  determining 
estimates  at  a  desired  epoch. 

1.3  Research  Approach 

The  basic  approach  to  this  research  took  many  steps.  First,  attempts  were 
made  to  learn  the  inputs  and  output  of  the  MSOFE  software.  This  program  was 
actually  the  Kalman  filter  "shell"  into  which  the  proposed  system  designs  had  to 
be  integrated.  The  filter  design  procedure  began  with  the  development  of 
mathematical  models  to  describe  the  actual  physical  system,  the  measurement 
system,  the  disturbance  process,  and  the  measurement  error  process  (Musick, 
1978).  In  the  strict  context  of  the  MSOFE  program,  the  program's 
documentation  distinguishes  between  the  "truth"  and  the  "system"  (Carlson  and 
Musick,  1990)  as  follows: 

Truth:  Defined  as  all  real-world  states  related  to  the  filter 
estimation  process.  This  model  is  the  starting  point,  from  which 
less  significant  states  are  systematically  eliminated  in  an  effort  to 
implement  a  more  (computationally)  practical  system. 


System:  A  subset  of  the  "truth,"  defined  by  all  significant  real- 
world  states  affecting  the  filter  estimation  process.  Included  are  the 
actual  values  of  the  filter  states  themselves  plus  higher-order 
measurement  error  states  that  significantly  affect  the  filter  states. 

Within  this  software,  the  "truth"  is  deemed  the  full-order  representation  of 
the  real  world,  the  "filter"  is  the  first-order  representation,  and  the  "system"  is 
the  second-order  representation3.  The  program,  MSOFE,  requires  adequate 
mathematical  representation  of  both  the  "system"  (also  called  the  "system  truth") 
and  the  "filter"  models.  The  main  focus  in  this  thesis  was  to  examine  the  effects 
on  the  "filter"  performance  due  to  higher-order  error  sources  present  in  the 
"system  truth"  model  (absent  in  the  simplified  "filter"  model).  The  challenge 
was  to  develop  a  "filter"  model  that  contained  a  sufficient  number  of  states  and 
sufficient  measurement  information  to  represent  adequately  the  physical 
phenomenon  of  interest,  here  the  GPS  process  and  measurement  models.  Thus 
the  model  used  for  filtering  was  a  much  simpler  model/subset  of  the  system  truth. 
Four  station-dependent  bias  sources  (the  station  coordinates  x,  y,  and  z;  and  the 
station  clock)  were  maintained  within  the  "filter"  model.  An  addition  of 
ambiguity  unknowns  was  added  when  satellite  phase  data  were  also  considered. 
The  "system,"  which  represented  the  real  world,  contained  combinations  of  the 
three  different  bias  types  listed  above. 

3Here,  the  "real  world"  refers  to  the  world  as  is  seen  from  day  to  day-a  state  of  being  which  is 
impossible  to  accurately  model.  Subsequently,  the  term  "real  world"  will  be  used  to  define  the 
model  corrupted  by  biases  (versus  the  "ideal  world"  with  no  biases).  The  distinction  will  be  made 
again  later  for  clarity. 


» 
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The  first  step  in  writing  software  code  to  utilize  the  MSOFE  program,  once 
the  details  of  the  program  were  deciphered,  was  to  develop  the  system  model  and 
the  correct  measurement  model  for  the  GPS  pseudorange  and  phase  range 
measurements.  Recall  that  both  the  "system"  and  the  "filter"  models  had  to  be 
distinguished  between  for  MSOFE  to  appropriately  estimate  the  desired 
quantities. 


» 


> 


The  incorporation  of  biases  into  the  pseudorange  and  phase  range 
measurement  models  was  a  step-wise  procedure.  The  "filter’s"  state  vector— the 
vector  of  estimated  quantities— contained  the  four  error  sources  inherent  in  any 
application  of  GPS:  the  station-dependent  biases,  listed  above.  Augmenting 
these  basic  four,  biases  were  added  to  the  "system  truth”  state  vector  to  better 
represent  the  GPS  case.  First  added  were  the  observation-dependent  biases  and 
integer  ambiguities. 


> 


> 


» 


Next,  tropospheric  delay  models  were  added  for  phase  and  range 
measurements.  The  model  for  calculating  the  delay  was  taken  from  a  currently 
used  Hopfield  model  (Goad  and  Goodman,  1974).  The  delay  caused  by  the 
ionosphere  was  next.  Values  for  this  bias  source  came  from  the  model  defined 
by  the  GPS  receiver  Interface  Control  Document  (ICD),  the  model  taking  its  inputs 
from  the  information  within  the  broadcast  ephemeris  (Leick,  1990). 


In  an  attempt  to  examine  different  measurement  weightings,  next  was 
created  a  measurement  noise  model  which  depended  on  the  elevation  of  the 


l 
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satellite  with  respect  to  the  station.  Two  different  models  were  used:  first,  with 
the  measurement  noise  dependent  on  the  tropospheric  delay  of  each  satellite,  and 
second,  on  the  ionospheric  delay  of  each  satellite4.  The  models  weighted  each 
error  source  based  on  the  satellite's  elevation,  the  lower  satellites  contributing 
less  weight  than  the  higher  ones.  Based  on  the  mathematical  model,  the  four 
station-dependent  unknowns  required  the  presence  of  [at  least]  four  satellites  to 
estimate  uniquely  their  values.  However,  only  the  addition  of  redundant  satellites 
would  affect  the  outcome  using  this  elevation-dependent  model.  Results  were 
expected  to  change  only  after  the  addition  of  a  fifth  satellite,  because  the 
weightings  for  the  four-satellite  case  have  no  influence  due  to  the  lack  of 
redundancy. 

The  next  iteration  explored  the  use  of  pseudorange  and  phase  range 
measurements  using  an  equal  measurement  uncertainty  of  ±  1  meter. 

Next,  the  satellite-dependent  biases  were  added  to  the  model. 
Corresponding  states  in  the  state  vector  were  (per  satellite)  two  trigonometric 
coefficients  and  one  satellite  clock  error  term. 

Finally,  the  issue  of  cycle  slips  was  investigated  by  introducing  a  simulated 
cycle  slip  into  the  set  of  "filter"  equations. 


4Both  tropospheric  delay  and  ionospheric  delay  are  functions  of  the  satellite's  elevation.  These 
models  will  be  called  "elevation-dependent"  models. 
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The  measurement  data  and  biases  in  this  thesis  were  all  simulated:  each  piece 
of  information  was  calculated  at  each  epoch.  Attempts  were  made  to  derive 
measurement  values  that  emulated  actual  observations  as  closely  as  possible.  Input 
for  calculating  the  data  came  from  actual  GPS  broadcast  ephemeris  information  and 
common  error  model  schemes.  The  actual  time  tags  "seen"  by  the  broadcast 
ephemeris  in  generating  orbit  information  were  incremented  in  order  to  simulate  the 
movement  of  the  satellites  during  the  test  time.  After  every  new  iteration,  the 
outputs  of  the  filter  were  examined.  Interest  lay  in  the  effects  of  these  iterations  on 
the  position  uncertainty  of  the  station/receiver.  Therefore,  the  results  were 
presented  in  the  form  of  graphs  showing  the  receiver's  position  uncertainty  versus 
time,  where  the  uncertainty  was  presented  as  an  unsealed  dilution  of  precision 
(DOP)  in  meters5.  The  results  from  each  iteration  were  compared  to  the  ideal  cases 
and  similar  combinations  of  biases.  Details,  results,  and  further  information  are 
presented  in  Chapter  IV,  Research  and  Conclusions. 


ft 
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5  Actually,  this  uncertainty  is  the  square  root  of  the  trace  of  the  filter's  error  estimation 
covariance  matrix  (covariance  of  estimated-minus-true  states).  DOP’s,  as  such,  are  in  units  of 

"meters  per  meter, "  and  are  defined  as  the  square  root  of  the  trace  of  the  filter's  error  estimation  ft 

covariance  divided  by  the  uncertainty  of  the  measured  variable.  Thus,  here,  the  uncertainty  was 
presented  in  meters. 
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CHAPTER  II 

THE  EXTENDED  KALMAN  FILTER 

2.1  Overview 

This  chapter  is  designed  to  provide  a  review  of  Kalman  filtering  definitions, 
techniques,  and  applications  to  this  thesis  topic.  Some  fundamental  knowledge  (by 
the  reader)  of  Kalman  filters  is  assumed.  The  Kalman  filter  theory,  as  it  pertains  to 
the  software  package  employed  in  this  thesis  research,  will  also  be  addressed.  If 
further  details  about  the  Kalman  filter  are  desired,  the  reader  may  obtain  helpful 
information  in  the  textbooks  by  Brown  (1983),  Gelb  (1992),  and  Maybeck  (1979). 

2.2  Introduction  to  the  Kalman  Filter 

According  to  Maybeck  (1979),  a  Kalman  filter  is  an  "optimal  recursive  data 
processing  algorithm"1  .  The  definition  may  be  broken  up  into  parts  to  help  clarify 
it. 

1The  concise  definition  of  a  Kalman  filter  differs  from  book  to  book.  Fundamentally,  the  Kalman 
filter  must  be  thought  of  as  a  concept  and  not  as  a  tangible  item  such  as  an  algorithm.  Maybeck's 
definition  expounds  upon  this  idea  by  implying  that  the  algorithm  is  the  means  by  which  the 
concept  is  carried  out. 


a.  "Optimal"  implies  that  the  filter  can  incorporate  all  information 
provided  to  it.  It  processes  all  available  measurements  to 
estimate  the  variables  of  interest  to  the  developer.  "Optimal" 
also  means  that  errors  are  minimized,  in  some  respect. 
According  to  Maybeck,  the  filter  is  able  to  utilize  knowledge  of 
system  and  measurement  dynamics,  statistical  descriptions  of 
any  noises  and  biases,  and  initial  conditions  of  the  variables  of 
interest  (ibid). 

b.  "Recursive"  implies  that  the  filter  algorithm  does  not  require  the 
storage  of  huge  amounts  of  data.  As  the  data  are  processed, 
the  updates  are  made,  and  then  the  processed  data  are  not 
needed  for  any  subsequent  updating.  Information  pertaining  to 
the  past  data  is  maintained  sequentially  within  the  update 
process2. 

c.  "Data  processing  algorithm"  implies  that  the  filtering  procedure 
is  realized  through  the  incorporation  of  discrete-time 
measurement  data,  rather  than  continuous  time  inputs.  The 
Kalman  filter  not-so-aptly  shares  its  name  with  the  type  of 
"black  box  filter"  known  to  the  electrical  engineer,  but  here,  the 
filtering  process  is  generated  more  typically  by  a  computer 

Recursive  adjustment  of  data  is  also  called  ’sequential'  adjustment,  in  which  measurements  are 
incorporated  sequentially,  as  they  become  available.  This  method  is  in  contrast  to  the  'batch* 
adjustment,  in  which  all  available  measurements  are  used  together,  processed  all  at  the  same  time. 
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program.  However,  the  main  idea  behind  the  two  types  of 
filters  is  the  same-imperfect  inputs  are  processed  into  outputs 
which  are  optimized,  with  respect  to  the  model  assumptions. 

Figure  2  shows  the  typical  application  for  a  Kalman  filter  (Maybeck,  1979).  A 
system  is  driven  by  control  variables,  which  are  affected  by  errors  and/or  biases. 
The  difficulty  in  "solving"  the  problem  at  hand  is  that  knowledge  of  the  system 
inputs  and  outputs  may  be  the  only  information  the  user  has  to  determine 


SYSTEM 
•  STATE 
ESTIMATE 


FIGURE  2 


TYPICAL  APPLICATION  OF  A  KALMAN  FILTER 
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the  state  of  the  system.  The  variables  of  interest  will  be  called  the  system  "states" 
and  may  or  may  not  be  directly  measured.  The  measuring  devices  are  the  means  for 
using  collectible  data  to  solve  the  problem  of  finding  the  estimates  of  the  variables. 

I 

Several  different  measuring  devices  may  be  available,  as  well.  However,  these 

measurement  devices  are  also  affected  by  error  sources,  usually  in  the  form  of  biases 

and  noises.  The  observations,  together  with  the  knowledge  (current  and  prior)  of 

the  system  and  measurement  devices,  are  combined  in  the  Kalman  filter.  Prior 

knowledge  may  include  the  a  priori  knowledge  of  the  system  description  and  the 

measurement  models.  The  output— the  desired  results— are  optimal  estimates  of  the  ( 

variables  of  interest,  more  accurate  than  the  estimates  based  on  the  individual 

measurements  (Stacey,  19°  1). 

► 

The  filtering  process,  as  was  indicated  above,  attempts  to  derive  the  optimal 
estimates  from  data  obtained  from  the  real  world-a  typically  noisy  and  error-laden 
environment.  The  filter  is  designed  to  propagate  the  conditional  probability  density  } 

of  the  desired  quantities,  meaning  the  shape  and  location  of  the  function  depend  on 
the  measurements  taken  (Maybeck,  1979). 

» 

A  further  assumption  with  respect  to  the  Kalman  filter  is  that  the  system  be 
described  by  a  linear  model,  and  the  system  and  measurement  noises  are  white  and 
Gaussian  (ibid)3.  The  "linear"  requirement  may  be  expanded  to  a  non-linear  system,  ^ 

provided  the  set  of  system  equations  is  linearized  about  some  nominal  point  or  value 

3  These  two  requirements  are  necessary  for  the  maximum  likelihood  interpretation  of  the  optimal 

estimate  to  hold  true.  I 


I 


not  too  far  away  (relatively)  from  the  truth.  This  linearized  model  then  becomes  an 
error  model,  and  filter  process  outputs  will  yield  error  estimations. 

The  "whiteness"  of  the  noise  means  that  the  noise  is  not  correlated  in  time- 
knowledge  of  the  noise  at  any  time,  t ,  will  not  allow  the  user  to  predict  the  noise  at 
time,  t+1.  The  actual  concept  of  the  "whiteness"  of  noise  is  just  an  idealized  means 
for  visualizing  the  noise's  "wideband"  behavior  in  a  "narrowband"  system. 
"Gaussian,"  on  the  other  hand,  pertains  to  the  amplitude  of  the  probability  density, 
taking  on  the  shape  of  a  bell-shaped  curve,  called  the  normal  probability  density 
function.  The  first-  and  second-order  statistics  of  a  noise  process— namely  the  mean 
and  variance— define  the  Gaussian  density. 

2.3  The  Kalman  Filter  Models 

Gelb  (1992)  states  that  before  the  filtering  process  can  be  successfully  applied, 
proper  modeling  and  realistic  performance  projections  must  be  addressed.  Three 
steps  must  be  incorporated  into  the  practical  application  of  estimation  theory: 

a.  design  and  computer  evaluation  of  the  "optimal"  system 
behavior; 

b.  the  design  of  suitable  "suboptimal"  system  with  cost  constraints, 
sensitivity  characteristics,  computational  requirements,  meas¬ 
urement  schedules,  etc.;  and 
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c.  the  construction  of  and  test  of  a  prototype  system,  making  final 
adjustments  or  changes  as  warranted. 

2.3.1  Basic  Statistical  Definitions 

The  expectation  of  a  random  variable,  x,  is  defined  as  the  sum  of  all  values  the 
random  variables  may  take,  where  each  value  is  weighted  by  the  probability  density 
with  which  the  value  is  taken  (ibid).  Also  called  the  mean  value  of  x,  the 
expectation  of  x  is  given  by 


E[x]=/::xf(x)dx  (2-1) 

where  f(x)  is  the  probability  density  function.  In  a  probabilistic  sense,  the 
expectation  is  the  value  to  which  the  average  of  a  number  of  observations  of  x  will 
tend  to  become  as  the  number  of  observations  increases. 

The  mean  squared  value ,  describing  the  distribution  of  x,  is  the  expectation  of 
x2,  written  as 


E[x2]=J::  x2f(x)dx  (2-2) 

and  the  root-mean-square  (rms)  is  the  square  root  of  the  estimated  E[x2].  The 
variance  (a2)  of  x  is  defined  as  the  expectation  of  the  deviation  of  the  random 
variable  from  its  mean  (ibid): 


» 
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oJ=J::(x-E[x])Jf(x)dx  (2-3) 

=  E[x2]-E[x]2 

The  square  root  of  the  variance  is  called  the  standard  deviation  of  the  random 
variable.  One  should  note  that  the  rms  value  and  the  standard  deviation  are 
comparable  only  for  the  zero-mean  random  variable. 


The  covariance ,  the  measure  of  how  the  error  in  one  random  variable  is 

I 

related  to  the  error  in  another,  is  the  expectation  of  the  product  of  the  joint 
deviations  from  their  respective  means,  written  as: 

E[(x  -  E[x])(y  -  E[y]>]  =  iZ  \Z  (x  -  E[x])(y  -  E[y])  f2(x,y)  dx  dy  1 

=  E[xy]  -  E[x]E[y]  (2-4) 

where  f2(x,y)  is  the  joint  probability  density  function  of  two  random  variables,  x 
and  y. 

> 

The  normal  probability  density  function,  which  characterizes  the  "Gaussian" 
process,  is  expressed  as 

f(x)=^exp[-^]  (2-5) 


I 


» 


» 
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where  p.  is  the  mean  and  o  is  the  standard  deviation  of  the  random  variable  x.  One 
may  write  the  n-dimensional  Gaussian  vector  s  as  in  Equation  2-64. 

S  ~  N  (Ul  2)  (2-6) 

Equation  2-6  represents  the  standard  notation  for  the  random,  Gaussian  (normal) 
vector,  x,  with  mean  ja  and  covariance  matrix  2. 

2.3.2  Kalman  Equations 

This  section  will  outline  the  equations  common  to  the  Kalman  filter 
mathematics  and  practical  application.  These  equations  are  employed  by  MSOFE. 

2.3.2.1  Transition  Matrix 

The  dynamics  of  a  continuous  system  may  be  represented  by  the  first-order 
differential  equation5 

x(t)  =  F(x(t),t)  x(t)  +  G(x(t),t)  w(t)  +  L(x(t),t)  n(t)  (2-7) 

where  a(t)  =  random  state  vector 

w(t)  =  random  forcing  function 


» 

> 


4The  boldface  notation  indicates  a  matrix,  while  the  underline  notation  indicates  a  vector. 

3 As  all  these  vectors  and  matrices  are  time-dependent,  the  reference  to  time  t  will  be  dropped  from 
subsequent  equations,  unless  needed  for  clarity. 


jj(t)  =  deterministic  (control)  input  vector 
F(x(t),t)  =  homogeneous  state  dynamics  matrix 
G(x(t),t)  =  system  process  noise  mapping  matrix 
L(x(t),t)  =  control  input  matrix. 

Typically,  the  relation  underlying  Equation  2-7  is  a  non-linear  ordinary  differential 
equation,  which  must  be  linearized  to  be  utilized  in  this  context.  This  equation  must 
be  known  together  with  stochastic  prior  information  used  in  the  linearization 
process,  if  Equation  2-7  is  derived  from  a  non-linear  equation,  in  order  for  the 
process  of  Kalman  filtering  to  hold  (more  about  the  linearization  of  this  equation 
later).  In  this  thesis,  the  GPS  problem  did  not  demand  any  controlling  inputs,  so  the 
deterministic  control  term  Lu  =  0,  and  Equation  2-7  for  this  case  becomes 

x  =  Fx  +  Gw  (2-8) 

Further,  the  mapping  matrix  G  is  typically  represented  as  the  identity  matrix  I,  as 
dictated  by  the  prior  information  used  for  linearization  of  the  differential  equation  (if 
necessary).  When  G  =  I,  the  system  noise  remains  unaltered  (unsealed)  as  white. 
Such  is  the  case  here.  Therefore,  Equation  2-8  is  written  as 

x  =  Fx  +  w  (2-9) 

The  differential  equation  given  above  must  be  sufficiently  general  such  that  any 
motion  can  be  described  through  one  of  its  solutions.  If  the  state  vector  and  forcing 
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functions  in  Equation  2-9  are  defined  at  time  t,  then  the  stales  at  any  other  time  can  9 

then  be  predicted. 

As  an  example,  one  may  substitute  into  Equation  2-9  a  particular  dynamics  * 

matrix  D  for  F  and  a  particular  white  noise  w'  for  w  to  get  Equation  2-10. 

x  =  Dx  +  w'  (2-10)  > 


Descritizing  Equation  2-10  and  then  through  rearranging,  one  obtains: 

-1*1  -  Dx  +  w’  (2-11) 

At 

xt+,  -xt  *  Dx  At  +  w’  At  (2-12) 

xt+1  “  +  Dx  At  +  w’  At  (2-13) 

-  (I  +  D  At )  x,  +  w'  At 

=>  Xfi  =  «&(t+l,t)  x,  +  u,  (2-14) 


> 


> 


where  u,  is  the  integrated  noise  vector  and  <X>(t+l,t)  represents  the  transition 
matrix.  The  transition  matrix  is  used  for  approximating  the  state  vector  at  time  t+1 
by  having  prior  knowledge  of  the  state  vector  at  time  t.  Equation  2-9  is  actually  the 


> 


> 
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equivalent  to  Equation  2-14;  the  former  the  continuous  representation,  and  the  latter 
the  discrete  representation6. 

2.3.2.2  Measurements 

Gelb  (1992)  equates  the  importance  of  measurements  to  the  observability  of 
the  system.  According  to  Gelb,  Gauss  defines  the  observable  system  as  one  which 
possesses  the  number  of  observations  that  are  absolutely  required  for  the 
determination  of  the  unknown  quantities.  The  basic  measurement  equation  may  be 
seen  in  Equation  2-15. 


z  =  Hx  +  v  (2-15) 

where  z  =  z(t)  =  the  measurement  vector 

H  =  H(x(t),t)  =  matrix  that  linearly  relates  x  to  z  (observation  matrix) 
v  =  white  measurement  noise  vector. 


^The  distinction  between  continuous  and  discrete  must  be  made.  Although  the  system  process  is  a 
continuous  one,  the  sampling  of  information  (measurements,  bias  determinations,  etc.)  is 
conducted  at  specified,  discrete  time  intervals  (Brown,  19S3).  In  this  thesis,  the  sampling  times 
are  intentional  and  under  the  control  of  the  programmer.  As  the  system  is  considered  continuous, 
the  Kalman  filter  equations  will  be  modeled  for  a  continuous-time  linear  system  with  discrete 
measurements. 
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2.3.2.3  Errors 


» 


ft 


*■ 


Recall  Equation  2-9,  assumed  now  to  be  a  linear  differential  equation  having 
been  written  in  terms  of  the  "whole  value"  state  variables.  However,  the  models 
employed  in  this  thesis  are  non-linear  and  must  be  linearized,  forming  error  models. 
Thus,  the  state  vector  written  in  terms  of  the  error  state  variables,  becomes 


6x  =  F  6  x  +  w 


(2-16) 

ft 


where 


Likewise,  for  the  "whole  value"  measurement  model.  Equation  2 


15,  the  measurement  error  model  becomes 


6  z  =  H  a  x  +  v  (2-17) 

ft 

where  Hh  -  . 

'  5xi 

ft 

In  Equation  2-18,  one  sees  the  representation  of  the  error  in  the  estimate 
vector  (with  the  "tilde"  -)  in  terms  of  the  estimated  (with  the  "hat"  a)  and  the 
actual  state  value. 


-  x  -  x 


(2-18) 

ft 


ft 
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The  covariance  matrix  Pof &t  will  therefore  be 

P-E[e,|J]  (2-19) 

assuming  that  x  is  an  unbiased  estimate  of  £.  To  cite  a  small  example  from  Gelb 
(1992),  given  the  error  estimate  as  the  two-state  vector  in  Equation  2-20, 

ix  ”  [  - 1  ]  (2*20) 

e2 

the  covariance  P  will  be  (under  the  assumption  that  the  expectation  of  each 
respective  state  is  zero) 


E{[ 


C|C2 


E[ef  ]  E[e,e2] 

Efe;] 


e2  ^  E[CjC2] 


] 


(2-21) 


An  important  item  to  note  is  that  the  square  root  of  the  trace  of  P  is  the  rms  length 
of  the  uncertainty  of  the  estimated  state  vector  x . 


The  covariance  for  the  white  noise  forcing  function  Gw  is  defined  as  (Gelb, 
1992) 


E[(GlwtXGlwi)T]  -  G(QtG^6(t-x) 


(2-22) 
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where  Q  is  the  process  noise  strength  matrix  and  5  is  the  Dirac  delta  function. 
Information  for  Q  comes  from  the  prior  information  used  for  the  linearization  of 
Equation  2-8.  If  one  lets  the  process  noise  distribution  matrix  G  be  the  identity 
matrix.  Equation  2-22  becomes 


E[w,w?]  =  Qt8(t-T) 


(2-23  a) 


T  Q  for  t,  =  t 
'  E[w(t  )w(t|)T]  =  {  ^  ,  J 

v  j  /  j  '  o  for  t,  *tj 


(2-23b) 


which  is  now  the  covariance  matrix  of  the  white  noise  w.  One  can  also  write  that 
the  mean  of  this  white  noise  vector  is 


4*  =  E[w]  =  0 


(2-24) 


Likewise,  for  the  white  measurement  noise  vector  v,  one  can  write  that  the  mean  is 
(Maybeck,  1979),  assuming  proper  modeling  in  Equation  2-15, 


=  E[v]  =  0 


(2-25) 


and  the  covariance  is 


T  R  for  t,  =  t 
E[v(t,)v(t,)T]  =  {  „  J 

1  ~v  J  1  1  0  for  t,  ^  tj 


(2-26) 


One  should  also  note  that  no  correlations  occur  between  w  and  v,  as  seen  in 
Equation  2-27. 


E[w(ti)v(tj)]-0 


(2-27) 


2.3.2.4  Update  Equations 

The  process  of  the  Kalman  filter  is  designed  to  "propagate"  (by  numerical 
integration)  the  error  state  vector  and  its  covariance  from  the  instant  in  time 
immediately  following  the  most  recent  measurement  update  t*  to  the  instant  in  time 
immediately  preceding  the  next  measurement  update  tr.,  (Stacey,  1991)7.  The 
values  for  the  state  vector  and  covariance  matrix  at  t'+l  are  the  predicted  values 
from  which  the  next  updates  are  made  (using  the  measurements).  The  following 
differential  equations  used  in  this  integration  are  given  (without  proof)  in  Equation 
2-28  and  2-29  (Maybeck,  1979). 


(2-28) 
(2-29) 

at  some  time  ^  with  the  initial  conditions 


x  -  Fx 

p-fp  +  pft+gqgt 


7The  notation  "super  indicates  the  condition  prior  to  the  measurement  update  (the  predicted), 
and  the  "super  +"  indicates  the  condition  after  the  update. 
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x-x(C) 

P  -  P(t*)  (2-30) 


as  provided  by  the  measurement  cycle  at  time  ts . 

After  propagation,  the  estimates  of  a  and  P  are  updated,  meaning  that  state 
estimates  are  revised  based  on  new  information  (Stacey,  1991).  Recall  Equation  2- 
15  for  the  measurement  model.  One  can  write  the  updated  state  vector  as 

x*  -  x  +  KAz  (2-31) 

where  Az  is  the  measurement  residual  and  K  represents  the  time -varying  Kalman 
filter  gain  matrix,  an  integral  element  in  the  update  equations  (to  be  explained  later). 
The  measurement  residual  may  be  written  in  terms  of  the  state  prior  to  the 
measurement  update  (x',  thought  of  as  the  predicted  state  vector),  changing 
Equation  2-31  to  Equation  2-32: 

X*  -x"  +  K(z-Hx*)  (2-32) 

Substituting  in  the  measurement  z ,  written  in  terms  of  the  actual  state  at  the  time  t. 
Equation  2-32  becomes 


» 


5> 


» 


» 


» 


» 


> 


» 


x*  -  \  +K[(Hx  +  v)-Hx’] 


(2-33) 

» 
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By  letting  Equation  2-34  represent  the  error  of  the  state  based  on  the  prior  value  of 
the  state, 

8x~=x'-x  (2-34) 

then  Equation  2-33  may  be  written  as 

x+  =  x_  +  K( — H5x~  +  v)  (2-35) 

Subtracting  x  from  each  side  and  defining 

8x+=x+-x  (2-36) 

as  the  error  of  the  state  based  on  its  updated  value,  one  gets 

5x+  =  5x  +  K(-H5  x"  +  y)  (2-37) 

Rearranging  Equation  2-37  gives  the  equation  for  the  state  error  vector  update: 

5x+  =  (I  -  KH)5x“  +  Kv  (2-38) 

Then,  the  covariance  of  this  updated  state's  error  is 

P+  =  (I  -  KH)P~(I  -  KH)T  +  KRKt  (2-39) 


I 


I 


* 
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where  R  is  defined  in  equation  2-25  and  P"  is  the  predicted  covariance  matrix  prior 
to  the  measurement  incorporation. 

The  Kalman  gain  matrix  K  should  be  discussed.  This  gain  is  a  means  for 
weighting  new  information  as  it  becomes  available  to  the  system.  This  new  added 
information  can  be  thought  of  as  causing  the  difference  between  the  actual 
measurement  and  the  filter's  estimate  of  this  measurement  (Stacey,  1991).  The 
Kalman  gain  matrix  is  designed  to  improve  the  estimate  of  the  state  vector  by  using 
information  contributed  by  the  observations.  The  gain  matrix  uses  knowledge  of 
measurement  noise  statistics  and  filter-computed  covariances  from  previous  epochs 
to  make  the  updates.  Without  proof,  the  Kalman  gain  matrix  for  a  continuous 
system/discrete  observation  model  is  (Gelb,  1992) 

K  =  P~(H~)t[H'P"(H~)t  +  R]'1  (2-40) 

where  H~  -  H(x').  This  value  of  K  is  the  precise  value  that  minimizes  the  trace  of 
the  P*  covariance  matrix  (ibid).  Substituting  Equation  2-40  into  Equation  2-39,  one 
gets  the  optimized  value  for  the  updated  estimation  error  covariance  matrix,  seen  in 
Equation  2-41. 


P*  -(I-KH)P' 


(2-41) 


ft 


* 


ft 


5> 


» 


ft 


» 


> 


> 


> 


> 


27 

2.3.2.5  Updated  Equations  with  Real-World  Errors 

The  previous  section  described  the  update  procedure  when  standard 
measurement  models  (with  noise  only)  are  used.  However,  in  the  real  world,  the 
measurement  models  and  state  equations  should  include  biases  due  to  real-world 
problems8.  The  more  realistic  measurement  model  may  be  written  as 

z  =  H  x  -  My  +  v  (2-42) 

where  y  =  the  vector  of  additional  random  error  states  (here,  called  biases), 

M  =  corresponding  mapping  matrix,  and 
the  negative  sign  is  used  only  as  a  convention. 

The  update  of  the  state  vector  may  now  be  written,  as  in  Equation  2-3 1,  as 

x*  =  x~  +  K(Hx  -  My  +  v  -  Hx~ ) 

=  x'  +  K(-H5x~  -  My  +  v)  (2-43) 

Following  the  same  procedure  as  in  Section  2. 3. 2. 4,  one  can  write  the  state  error 
vector  update  equation  as 


5x*  =  (I  -  KH)5x~  -  KMy  +  Kv 


(2-44) 


8Recall  the  definition  of  "bias"  here  is  not  defined  in  the  statistical  sense,  but  implies  a  random 
influence  on  the  GPS  measurement  which  must  be  modeled  (systematically)  and  accounted  for. 


and  the  covariance  of  this  updated  state's  error  is 

P*  =(I-KH)P'(I-KH)T+KRKT+KMPTMTKT  (2-45) 

where  PT  is  the  covariance  of  the  error/bias  (real-world)  vector;  E[y]  =  0;  and  no 
correlations  exist  between  y  and  v  or  w,  respectively. 

One  should  note  a  very  important  point  in  Equation  2-45.  The  last  term 
represents  the  addition  to  the  state  vector's  uncertainty  due  to  the  presence  of  the 
unmodeled  errors  in  the  observation  equations.  MSOFE  enables  the  user  to  look  at 
both  of  the  following9: 


a.  The  output  of  the  filter,  i.e.,  the  filtering  process's  estimate  of 
the  error  it  commits  in  estimating  the  entire  filter  state  vector 
(represented  by  Equation  2-39);  and 

b  The  "true"  state  vector  estimation  error,  i.e.,  the  "true” 
covariance  of  error  in  the  filter-computed  state  as  a  random 
process  driven  by  the  true,  real-world  system  state  (represented 
by  Equation  2-45). 

Furthermore,  one  must  assume  here  that  the  y  vector  is  a  vector  of  "random 
biases"  which  affect  the  GPS  measurement  and  does  not  contain  quantities  to  be 


9The  distinction  here  is  one  strictly  imposed  by  MSOFE.  Recall  that  MSOFE  distinguishes 
between  the  "system  truth"  (or  equivalently  "system")  and  the  "filter"  state  vectors. 


estimated  by  the  filtering  process.  In  turn,  the  covariance  matrix,  PT,  will  be 

considered  bounded.  This  assumption  will  allow  one  to  consider  only  defined 
covariances  and  ignore  the  unique  situation  when  P?  becomes  unbounded  (i.e., 

approaches  infinity). 
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CHAPTER  III 

GLOBAL  POSITIONING  SYSTEM 

» 

3.1  Basic  Theory 


In  this  chapter,  some  of  the  basic  concepts  behind  the  Global  Positioning  System 
(GPS)  will  be  summarized,  beginning  with  a  short  background  through  to  a  summary 
of  the  measurement  model  employed  in  the  research  for  this  thesis.  A  moderate 
amount  of  information  will  be  assumed  already  known  by  the  reader.  The  topics 
represented  here  play  a  roll  in  the  subsequent  chapter  dealing  with  the  actual  research. 
For  further  details  about  GPS  theory,  the  reader  may  obtain  helpful  information  from 
the  textbooks  by  Leick  (1990)  or  Wells  (1987). 

3.1.1  Background 

The  concept  for  the  GPS  is  to  serve  the  United  States'  military  community.  The 
United  States  Department  of  Defense  (DOD)  designed  the  system  to  be  a  military 
navigation  system  based  on  a  world-wide  coverage  of  a  constellation  of  24  satellites, 
or  space  vehicles  (SV)  (Stacey,  1991).  The  design  for  the  satellites'  orbit  and  earth 
coverage  intended  to  provide  24-hour  weather,  navigation,  timing,  and  surveying 
capabilities  (Leick,  1990). 


> 
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The  designed  constellation  of  satellites  was  originally  to  give  the  best  possible 
coverage  over  the  largest  area  of  concern  (in  the  military  field,  the  two  areas  of 
concern  were  the  southwestern  continental  United  States  and  off  the  northeast  coast 
of  the  United  States  in  the  North  Atlantic  Ocean).  The  satellites  are  positioned  in 
orbits  approximately  20,000  kilometers  above  the  surface  of  the  earth.  The  satellites 
have  a  orbital  period  around  the  earth  of  12  sidereal  hours.  Due  to  the  differences 
between  the  sidereal  and  solar  time,  the  satellites  appear  almost  4  minutes  earlier  each 
day,  the  coverage  advancing  by  approximately  30  minutes  each  week  (Wells,  1987). 

The  average  designed  life  span  of  each  satellite  ranges  from  5  to  7.5  years.  The 
satellites  are  equipped  with  solar  panels  which  charge  batteries  during  periods  of 
sunlight  Some  basic  functions  of  the  satellites  are  (ibid)  . 

(a)  To  receive  and  store  information  transmitted  by  the  system  operators; 

(b)  To  conduct  limited  data  processing  using  an  on-board  microprocessor; 

(c)  To  maintain  accurate  time  by  means  of  several  highly  accurate  oscillators; 

(d)  To  .ransmit  information  to  the  user  through  propagated  signals; 

(e)  To  maneuver  using  thrusters  controlled  by  the  system  operators. 

Due  the  height  of  the  orbit  of  the  satellites,  approximately  three  times  the  radius 
of  the  earth,  any  particular  satellite  may  be  visible  to  the  user  for  approximately  5 
hours  each  revolution.  Based  on  the  configuration,  time  of  day,  and  location  on  the 
earth,  the  user  may  "see"  4  to  10  satellites  above  the  horizon  at  any  one  time  (ibid). 
However,  due  the  geometry  of  the  satellites  with  respect  to  the  station  (user),  the 
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ability  to  recover  positions  will  vary  from  time  to  time.  In  fact,  the  geometry  of  the 
satellites  with  respect  to  a  receiver  plays  a  major  role  in  the  user's  ability  to  recover 
station  position  accurately.  Think  of  the  satellites  in  the  sky  and  the  receiver  on  the 
ground  forming  a  polyhedron;  the  larger  the  volume  of  the  polyhedron,  the  smaller 
the  recovered  position  uncertainty  (the  ideal  geometry  or  positioning  of  the  satellites 
over  the  receiver  places  one  satellite  overhead  and  the  other  three  on  the  horizon,  120 
°  apart  in  azimuth  This  geometry  (ideally)  allows  for  the  best  recovery  of  the 
receiver's  position  (ibid)).  As  the  satellites  move  through  the  sky,  this  geometry  will 
change,  as  will  the  volume  of  the  "polyhedron"  and  the  recovered  position 
uncertainty 

3.1.2  Signal  Structure 

The  fundamental  frequency  for  transmission  of  information  by  the  satellites  is 
defined  as  f0  =  10.23  MHz.  Each  satellite  transmits  two  carriers  at  two  different 

frequencies: 


L,  =  154  fo  =  1575.42  MHz  (3-1) 

L:  =  120  f„  =  1227.60  MHz  (3-2) 
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The  L,  is  modulated  with  two  types  of  code  (P-code  and  C/A-code),  and  the  L,  is 
modulated  with  only  one  type  of  code  (P-code);  both  contain  the  navigation  message. 

The  codes  are  called  pseudo-random  noise  (PRN)  codes  and  consist  of  a  series 
of  +  l's  and  -l's  emitted  at  particular  frequencies  If  the  code  value  is  -1,  the  carrier 
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phase  is  shifted  by  180°;  if  the  code  is  +1,  the  carrier  is  unchanged.  The  phase  shift  of 
the  signal  may  be  seen  in  the  following  picture,  indicated  by  the  arrow  in  Figure  3. 
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FIGURE  3 

REPRESENTATION  OF  PHASE  SHIFT  OF  CARRIER  PHASE 


The  C/A-code  (or  coarse/acquisition  code,  also  more  recently  called  Standard 
Positioning  Service  (SPS))  is  the  code  available  to  the  civilian  community.  The  C/A- 
code  is  a  relatively  short  code  of  1023  bits,  with  a  duration  of  I  millisecond  (msec) 
and  an  emission  (chipping)  rate  of  1.023  Mbps  (megabits  per  second)  (Spilker,  1980). 
The  chipping  rate  may  be  thought  of  as  the  rate  at  which  the  phase  is  shifted,  the 
means  for  determining  the  beginning  of  a  piece  of  data  (e  g.,  the  navigation  message). 
Each  satellite  transmits  mutually  exclusive  C/A  codes,  making  it  possible  to 
distinguish  signals  received  simultaneously  from  the  satellites  (Leick,  1990).  The 
C/A-code  is  only  available  on  the  L,  carrier. 

Considering  the  speed  of  light,  c*3x  10®  m/s,  the  length  of  the  C/A  signal  is 
approximately  300  meters.  One  can  expect  the  uncertainty  in  position  associated  with 
this  code  to  be  approximately  1%  of  the  code  length,  or  *  3  meters  (ibid). 


The  P-code  (or  precise  code,  also  more  recently  called  Precise  Positioning 
Service  (PPS))  emits  signals  at  10.23  Mbps,  repeating  itself  every  37  weeks  (267 
days).  This  code  is  initialized  at  the  beginning  of  each  week  (the  "GPS  week"  is 
defined  as  the  time  from  midnight  Saturday  to  Sunday  (ibid)),  but  is  not  mutualy 
exclusive  as  is  C/A  code.  However,  each  satellite  is  assigned  its  own  weekly  portion 
of  the  code,  and  some  weeks  will  not  be  utilized,  because  there  are  less  than  37 
satellites  in  the  GPS.  The  P-code  is  available  on  both  the  L,  and  the  L2  carriers. 
With  an  approximate  code  length  of  30  meters,  the  uncertainty  of  the  position 
recovery  is  a  0.3  meters  =  30  centimeters. 

A  third  type  of  code,  the  Y-code,  is  similar  to  the  P-code,  but  is  classified  by  the 
DOD  and  is  usable  only  to  those  users  who  have  access  to  the  encryption  key.  The 
intent  of  this  type  of  code  is  to  provide  a  secure  code  to  which  other  users  cannot 
deny  its  access.  This  method  of  selectively  providing  access  to  users  is  called  cinti- 
spoofmg  (AS)1. 

Both  carrier  frequencies  transmit  the  navigation  message.  This  message  is  a  50- 
Hz  stream  of  data  designed  to  notify  the  user  of  the  position  and  other  pertinent 
information  about  the  satellite.  The  message  contains  information  on  the  ephemerides 
of  the  satellites,  the  GPS  time,  the  behavior  of  the  clock,  and  any  system  status 
messages. 


1  This  type  of  code  will  not  be  addressed  further  in  this  thesis,  but  is  mentioned  only  as  background 
information. 


The  DOD  has  chosen  to  maintain  controls  on  the  users  of  the  GPS  by  a  method 
called  selective  availability  (SA).  When  desired,  the  Master  Control  Station  (the 
central  location  for  updating  broadcast  and  satellite  orbit  information)  located  in 
Colorado  Springs,  Colorado,  will  transmit  degraded  information.  The  broadcast 
information  may  include  "orbit  perturbations"  which  will  lead  to  degraded 
ephemerides.  The  atomic  clock  times  may  also  be  perturbed  (dithered)2. 

3.1.3  Pseudorange 

The  first  and  most  common  type  of  measurement  which  can  be  made  by  the 
GPS  receiver  is  the  pseudorange.  In  the  simplest  terms,  the  range  is  the  distance 
between  the  satellite  and  station  (receiver).  As  this  measurement  is  impossible  to 
explicitly  measure  with  ordinary  measuring  devices  (due  to  the  motion  of  the  satellites 
and  the  possible  motion  of  the  receiver),  the  GPS  receiver  uses  the  timing  of  signals 
to  arrive  at  the  satellite-station  distance.  The  time  the  signal  takes  to  travel  from 
satellite  to  station  is  given  by 

At  =  tR-ts,  (3-3) 

where  the  "R"  and  "S"  denote  "receiver"  and  "satellite,"  respectively.  Its  equivalent 
"distance"  would  be  cAt,  where  c  is  the  speed  of  light.  This  value  would  be  the  exact 
range  if  neither  the  clock  in  the  satellite  nor  the  clock  in  the  receiver  were  in  error. 
Alas,  this  exactness  does  not  exist  in  the  real  world.  We  expect  both  clocks  to  be  in 

2Sclcclivc  availability  serves  to  degrade  instantaneous,  real-time  position  calculations,  but  post¬ 
processing  methods  could  allow  a  more  accurate  determination  of  positions. 
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error  by  some  amount,  so  neither  clock  contains  the  correct  time.  The  measured 
time,  t,  on  either  the  satellite  clock  or  receiver  clock  can  be  given  by 

t  =  t  *  +8t  (3-4) 

where  t*  =  the  correct  (true)  time,  without  errors,  and 

St  =  the  correction  necessary  due  to  the  inaccuracies  in  the  clocks. 

Therefore,  the  actual  range  measured  by  the  GPS  receiver,  barring  any  other 
error  sources  or  biases  and  using  the  definition  for  range  above,  will  be 

R  =  ctR  -cts  (3-5) 

=>  R  =  c(tR  +5tR)-c(ts*  +8ts) 

=  c(t*  -ts*)  +  c8tR  -c8ts 

=  p  +  c5tR  -  c5ts  (3-6) 


with  p  representing  the  true  geometric  distance  between  satellite  and  receiver.  The 
value,  R,  is  called  the  pseudorange  because  of  the  presence  of  clock  errors  and  their 
addition  into  the  calculation  of  the  range.  To  be  more  precise,  Wells  (1987)  defines 
the  pseudorange  as  the 

time  shift  required  to  line  up  a  replica  of  the  code  generated  in  the 
receiver  with  the  received  code  from  the  satellite  multiplied  by  the 
speed  of  light 
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where  the  "time  shift"  is  the  difference  between  the  time  of  the  received  signal  and  the 
emitted  signal,  measured  in  the  receiver's  and  the  satellite’s  time  frames,  respectively. 
The  navigation  message  from  each  satellite  contains  the  necessary  information  for  the 
evaluation  of  the  satellite  clock  offset,  8ts,  in  Equation  3-6. 

The  distance  could  also  be  found  by  calculating  the  length  of  the  vector  from 
satellite  to  receiver  in  Equation  3-7, 

Pr  =yl(xs-*R)  +  (ys-yR)  +  (zs-zR),  (3-7) 

where  pR  =  range  from  satellite,  S,  at  epoch  ts  to  receiver,  R  at  epoch  tR;  and 

x,  y,  z  (normally)  denote  the  earth-centered,  earth-fixed  coordinates  of  the 
satellite  and  receiver.  Therefore,  including  clock  errors,  the  pseudorange  becomes: 

R|  =pR +c5tR -c5ts +eR  (3-8) 

where  p|  is  found  using  Equation  3-7, 

c5tR  -c5ts  is  the  clock  correction  term  from  Equation  3-6,  and 
eR  is  the  noise  term. 

3.1.4  Phase  Range 

The  second  type  of  measurement  we  can  extract  from  a  GPS  receiver  is  the 
phase.  According  to  Wells  (1987),  the  carrier  beat  phase  is  "the  phase  of  the  signal 
which  remains  when  the  incoming  Doppler-shifted  carrier  is  differenced  (beat)  with 


38 


the  constant  frequency  generated  in  the  receiver."  In  simple  terms,  the  phase  is  the 
number  of  cycles  the  carrier  travels  from  satellite  to  receiver3.  The  phase  range  is 
therefore  the  phase  converted  to  an  equivalent  distance  in  meters  (by  multiplying  by 
the  wavelength  of  the  signal).  Here,  phase  range  measurements  will  be  used.  Just  as 
with  the  pseudorange,  the  phase  range  is  impossible  to  measure  explicitly.  The 
equation  for  the  phase  range  measured  by  the  receiver  is 

(Pr(^r)  =  <P  0t)  —  (Pr(Jr)  (3*9) 

where  tR  =  the  time  at  the  receiver,  and 

tT  =  the  time  at  the  transmitter  =  tR  -  where  p  =  the  distance  traveled  in 
a  vacuum  and  c  =  the  speed  of  light. 

This  equation  for  phase  ranges  differs  from  that  of  pseudoranges  because  the 
latter  uses  the  differences  between  times,  and  the  former  differences  phase  range 
values  at  the  satellite  and  receiver.  The  phase  is  made  up  of  a  whole  number  of 
wavelengths  between  satellite  and  receiver  plus  the  fraction  of  wavelength  remaining. 
It  is  this  fraction  which  is  actually  measured  by  the  receiver;  the  whole  number 
(integer)  is  the  more  difficult  value  to  determine  Therefore,  one  may  think  of  the 
phase  range  as  (without  the  presence  of  any  other  biases) 


*PrOr)—  V  (^t)  VrOr)  Pr  (3-10) 


3Ho\vever,  this  is  not  what  is  measured  on  the  GPS  receiver.  This  will  be  explained  later. 
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where  the  first  two  terms  represent  the  non-synchronicity  of  the  two  clocks  expressed 
in  wavelengths,  and  p|  is  the  geometric  range  from  Equation  3-7.  If  the  clocks  were 
perfectly  aligned,  the  measured  phase  range  would  exactly  equal  the  geometric  range 
(barring  the  presence  of  any  other  biases  or  errors). 

Recall  the  values  for  the  frequencies  of  the  carrier  and  the  two  codes.  The 
carrier  wavelength  is  much  shorter  than  either  of  the  codes'  wavelengths.  Using  the 
same  logic  used  above  for  pseudoranges,  one  can  expect  the  precision  of  the  phase 
ranges  to  be  much  better  (higher)  than  that  of  the  pseudoranges.  For  the  L,  carrier, 
the  wavelength  is  approximately  19  cm.  Using  the  typical  value  of  1%,  the  phase 
measurements  can  be  made  to  approximately  1%  of  the  wavelength,  or  approximately 
0.2  cm  =  2  mm.  As  good  as  this  appears  on  the  surface  to  the  potential  users,  special 
problems  plague  the  phase  measurement  technique,  as  will  be  mentioned  a  bit  later. 

3.2  Measurement  Biases 

Because  the  real  world  is  continually  plagued  with  imperfections,  the 
measurements  attained  from  the  typical  GPS  receiver  will  not  be  representations  of 
the  exact  range  between  satellite  and  receiver.  Three  types  of  biases4  will  be  covered 
in  this  thesis: 


4Rcca!l  that  this  thesis  defines  the  "bias"  as  any  random  effect  on  the  GPS  measurement  which  must 
be  (systematically)  modeled  and  correct  for  in  order  to  obtain  an  accurate  range  measurement.  This 
definition  may  conflict  with  others'  notation  for  the  "bias"  as  a  non-random  quantity.  The  reader 
must  keep  these  differences  in  notation  in  mind. 
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a.  Station-dependent  biases, 

b.  Observation-dependent  biases, 

c.  Satellite-dependent  biases. 

3.2.1  Station-Dependent  Biases 

As  was  mentioned  above  in  Section  3.1.3,  the  receiver  clocks  contain  some 
error.  One  must  account  for  this  error  to  ensure  the  clock's  bias  does  not  affect  the 
range  values.  Similarly,  the  data  file  taken  from  the  receiver  gives  only  preliminary 
information  about  the  coordinates  of  the  station;  these  coordinates  are  just 
approximations.  In  fact,  these  four  values  (station  clock  and  the  three  station 
coordinates)  are  frequently  the  desired  quantities  to  be  extracted  from  a  series  of  data 

3.2.2  Observation-Dependent  Biases 

Three  observation-dependent  biases  exist  which  cause  much  consternation  in  the 
GPS  community:  tropospheric  delay,  ionospheric  delay,  and  carrier  phase  integer 
ambiguity  All  three  types  will  be  addressed  in  this  thesis.  The  first  and  second  types 
exist  because  of  the  presence  of  the  GPS  signal's  propagation  media. 
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3.2.2. 1  Propagation  Media 

Electromagnetic  radiation,  regardless  of  its  source,  is  affected  by  the 
propagation  medium.  Based  on  the  governing  laws  of  rays  through  a  medium,  the 
rays,  in  this  case  the  satellites'  signals,  are  said  to  refract  due  to  the  presence  of  the 
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medium.  This  bending  is  accompanied  by  a  delay  in  the  travel  time  of  the  rays.  The 
bending  action,  per  se,  does  not  effect  the  GPS  community  as  much  as  the  time  delay 
in  the  signals.  The  presence  of  the  media  physically  lengthens  the  true  path  of  the 
signal  with  respect  the  geometric  distance  from  the  satellite  to  the  receiver,  seen  in 
Figure  4  (Goad  and  Goodman,  1974).  The  lengthened  path  means  the  signal  has 
taken  a  longer  amount  of  time  to  reach  the  satellite  than  would  have  been  taken  in  a 
vacuum.  This  time  delay  biases  the  GPS  measurements  and  must  be  modeled  to 
remove,  or  at  least  suppress,  the  effects. 


ft 


ft 


ft 


ft 

ft 


ft 

FIGURE  4  * 

RANGE  VECTOR  VERSUS  TRUE  PATH 
Effect  of  Propagation  Media  on  GPS  Signals 

ft 
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The  propagation  media  for  a  typical  GPS  signal  are  the  troposphere  and  the 
ionosphere.  The  troposphere  is  the  portion  of  the  atmosphere  extending  from  the 


earth's  surface  to  about  40  km.  The  ionosphere  extends  from  approximately  50  km  to 
1000  km  in  altitude  (Wells,  1987).  The  two  media's  effects  on  the  GPS  signals  are 
quite  different,  as  are  the  means  for  modeling  them.  Both  will  be  discussed  here. 

3.2.2.1.1  Troposphere 

For  frequencies  below  30  MHz,  the  troposphere  is  said  to  behave  as  a 
nondispersive  medium.  A  nondispersive  medium  is  one  in  which  the  refraction  of  the 
signal  is  independent  of  the  frequency  of  the  signals  being  transmitted  through  it 
(Leick,  1990).  This  delay  is  very  dependent  on  the  pressure,  temperature,  and 
humidity  of  the  medium  (ibid).  Because  the  true  path  is  lengthened  (the  transmission 
time  of  the  signal  has  increased),  any  correction  due  to  the  bias  must  be  subtracted 
from  the  observed  GPS  measurement  (pseudorange  or  phase  range). 

The  troposphere,  as  in  any  atmospheric  medium,  is  made  up  of  a  dry  and  a  wet 
part.  Both  parts  are  modeled  to  predict  the  medium's  effect  on  the  signal  The  dry 
portion  of  the  troposphere  makes  up  about  90%  of  the  total  refractive  effect.  This 
portion  is  commonly  estimated  using  surface  pressure  and  temperature  data  and  is 
known  to  an  accuracy  of  about  2%  (Goad  and  Goodman,  1974). 

The  wet  portion,  although  only  10%  of  the  total  refractive  effect,  is  much  more 
difficult  to  model.  This  portion  is  dependent  on  the  conditions  along  the  transmission 
path  of  the  signal,  often  independent  of  the  surface  conditions.  The  model  must 
include  information  about  the  water  content  of  the  atmosphere,  the  temperature,  the 
altitude,  and  the  elevation  angle  of  the  signal-some  of  which  are  very  difficult  to 
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model.  The  major  problem  with  the  modeling  of  the  wet  part  lies  in  having  accurate 
knowledge  of  the  water  vapor  content.  Localized  water  "pockets"  in  the  atmosphere 
along  with  atmospheric  turbulence  create  more  problems  in  the  modeling  of  the  wet 
contribution.  Various  models  exist  which  model  this  component,  taking  all  the  above 
difficulties  into  account. 

The  overall  model  for  the  contribution  of  the  troposphere  to  the  GPS  signal  is  a 
combination  of  the  wet  and  dry  parts.  Based  on  the  information  found  in  three 
different  studies5,  the  average  correction  for  the  troposphere  based  on  the  satellite's 
elevation  with  respect  to  the  station  may  be  found  in  Table  1  (ibid). 
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TABLE  1 

AVERAGE  TROPOSPHERIC  CORRECTION 
Based  on  Elevation  Angle 
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Elevation  Angle  (decrees) 

Average  Correction  (meters) 

90 

2.40 

45 

3.39 

30 

4.78 

20 

6.95 

10 

13.34 

7.5 

17.34 

5 

24.52 

2 

45.61 

0 

88.28 

-  ! 

5TIic  three  studies  referred  to  arc  those  of  Saastamoincn,  Saastainoincn/Marini.  and  Hopficld 
Complete  references  may  be  found  in  Goad  and  Goodman  (1990). 
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3.2.2. 1.2  Ionosphere 


In  the  earth's  ionosphere,  the  sun  affects  a  fraction  of  the  gas  molecules  by 
ionizing  them  with  ultra-violet  radiation,  causing  the  release  of  free  electrons.  Just  as 
with  any  electromagnetic  signal  propagating  through  an  ionized  medium,  the  GPS 
signals  are  affected  by  the  non-linear  dispersive  characteristics  of  the  medium.  Unlike 
the  nondispersive  characteristic  of  the  troposphere,  the  ionosphere  is  frequency- 
dependent,  and  the  modulations  on  the  carrier  and  the  carrier  phases  are  affected 
differently. 

The  dispersive  effects  of  the  ionosphere  may  be  visualized  by  examining  the 
dispersion  curve  in  Figure  5  (Wells,  1987). 


FIGURE  5 

IONOSPHERIC  DISPERSION  EFFECTS  OF  THE  MEDIUM 
Ionospheric  Propagation  Wave  Number  vs  Angular  Frequency 
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At  any  point  in  the  curve,  the  slope  of  the  line  joining  the  origin  to  that  point  is  the 
phase  velocity ,  vt  =  where  co  =  27tf  is  the  angular  frequency  of  the  signal  in  the 

ionosphere,  and  p  =  2  7^/  js  the  ionosphere  propagation  wave  number.  The  group 
velocity ,  v§  =  is  found  by  the  local  tangent  slope  at  the  point.  The 

corresponding  indices  of  refraction  are 


n 


P 


and 


(3-11) 


where  c  =  speed  of  light  in  the  vacuum.  The  respective  indices  of  refraction  for  the 
ionosphere  are  (in  a  first-order  approximation) 


and 


(3-12) 


where  N  =  ionospheric  electron  density  in  electrons/m  ’  at  the  time  and  place  of 
evaluation  and  a  =  constant  (Wells,  1987).  The  change  in  path  length  due  to  the 
ionosphere  is  given  by 


Ap  =  |(n-l)ds  (3-13) 

which  represents  the  integration  along  the  propagation  path  The  GPS  signals  are 
dependent  on  the  group  index  of  refraction,  so  the  ionospheric  group  delay  is 
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aN 

Ap  =+ — — 

r  I  ^2 


(3 -14a) 


The  effect  on  the  phase  measurements  is  opposite  in  sign  to  the  range  effect  seen  in 
Equation  3- 14a;  the  ionospheric  phase  delay,  is,  therefore. 


aN 

Ap  =-- 
^*2 


(3- 14b) 


where  Nr  in  both  equations  is  the  total  electron  content  along  the  propagation  path  in 
electrons/m 

The  ionospheric  effect  also  sees  daily  variations,  the  minimum  between  midnight 
and  early  morning,  and  the  maximum  around  local  noon,  dependent  on  the  activity  of 
the  sun  (Leick,  1990).  Depending  on  these  factors,  the  total  ionospheric  effect  could 
range  from  centimeters  to  tens  of  meters  (Wells,  1987). 


To  reiterate,  the  most  important  effects  of  the  ionosphere  on  the  GPS  signal  are 
(Leick,  1990) 


a.  the  retardation  of  the  modulation  on  the  carrier  wave  (i.e.,  the 
ionospheric  time  delay),  causing  an  increase  in  the  apparent  path 
of  the  signal  (pseudorange);  and 


b  the  advancement  of  the  carrier  phase,  causing  a  decrease  in  the 
apparent  phase  of  the  signal  (phase  range). 


3.2.2.2  Carrier  Phase  Integer  Ambiguity 
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When  the  GPS  receiver  is  first  turned  on  to  take  measurements,  the  cycle 
counter  will  read  some  arbitrary  initial  integer  value.  If  the  satellite's  signal  is 
continually  received  by  the  station  during  the  measurement  period  (i.e.,  no  "loss  of 
lock"  occurs  to  disrupt  the  flow  of  information  from  satellite  to  receiver),  this  integer 
will  not  change.  Think  of  the  phase  range  measurement  Equation  3-10  as  now 
possessing  an  ambiguity  bias,  shown  in  Equation  3-15, 

» 

<PR(tR)  =  <f>S(tT)-<MtR)-pR  +  NX  (3'15) 


where  NX  is  the  integer  times  the  wavelength  of  the  carrier.  If  the  world  were  ideal, 
the  clocks  would  be  synchronous,  and  the  two  quantities  [rps(tT )  —  cpR  (tR  )]  and 
(the  reading  on  the  GPS  receiver)  would  equal  zero.  The  geometric  distance, 
therefore,  would  exactly  equal  the  integer  ambiguity,  NX  The  term,  NX,  represents 
the  integer  number  of  counts  between  the  satellite  clock  and  receiver  clock. 
However,  the  receiver  measures  the  fractional  portion  of  the  signal;  only  this 
fractional  part  is  meaningful  to  the  GPS  user. 
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Actually  determining  the  integer  ambiguity  value  can  be  very  difficult,  if  not 
impossible  It  is  difficult  to  locate  the  exact  cycle  of  the  carrier  whose  phase  is  being 

> 

measured  (Wells,  1974).  Although  not  meaningful,  this  integer  value  must  be 
resolved  to  properly  process  phase  data  One  way  to  eliminate  its  contribution  is  to 
difference  two  measurements  which  have  the  same  ambiguity  value,  possibly  at 

» 

consecutive  epochs,  canceling  out  the  ambiguity  from  the  problem,  altogether. 
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3.2.3  Satellite-Dependent  Biases 

Two  types  of  satellite-dependent  biases  exist:  orbit  representation  biases  and 
satellite  clock  biases.  Clock  errors  are  present  due  to  biases  in  the  models  for  the 
clocks  in  the  broadcast  message  (e  g.,  it  is  impossible  for  the  clocks  to  be  perfectly 
aligned  to  GPS  time).  Orbit  biases  exist  because  of  the  errors  inherent  in  the  satellite 
ephemeris  information  (e  g.,  the  satellite  is  not  where  the  broadcast  message  indicates 
that  it  is).  Errors  in  the  ephemeris  may  be  eliminated  by  access  to  better  orbit 
information  The  GPS  user  is  not  often  lucky  enough  to  receive  the  most  accurate 
information,  and  as  such,  must  compensate  for  these  biases.  Forces  not  under  the 
user's  control  continually  act  on  the  satellites,  rendering  what  information  the  user 
does  have  as  inaccurate.  These  forces  are  also  very  difficult  to  measure  from  any 
location  on  the  earth's  surface. 

Orbital  biases  may  naturally  be  expressed  in  the  form  of  radial,  along  track,  or 
out  of  plane,  as  seen  in  Figure  6.  Using  orbital  data  provided  by  the  broadcast 
ephemeris,  satellite  positions  can  expect  a  typical  accuracy  of  about  10-20  meters 
(Wells,  1987). 


RADIAL 
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FIGURE  6 
ORBITAL  BIASES 


Here,  the  orbital  bias  introduced  into  the  measurement  models  will  be  a 
combination  of  the  three  different  types.  The  combined  bias  may  be  thought  of  a 
periodic  oscillation  of  the  satellite  along  its  orbit.  These  oscillations  are  fractions  of 
the  sine  and  cosine  of  the  angular  position  of  the  satellite  as  a  function  of  time, 
defined  by  the  satellite's  Keplerian  elements,  as  in  Equation  3-16.  The  value  for  Ar 
indicates  the  change  in  the  radius  of  the  orbit  based  on  the  trigonometric  oscillations 
of  the  satellite. 


Ar  =  Ssin4>  +  Ccos({>  (3-16) 

where  <J>  =  argument  of  perigee  (to)  +  true  anomaly  (f) 

S  and  C  =  coefficients 
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The  argument  of  perigee  and  the  true  anomaly  place  the  satellite  in  the  orbital  plane 
at  a  particular  epoch,  as  seen  in  Figure  7. 


FIGURE  7 

PLACEMENT  OF  THE  SATELLITE  IN  THE  ORBITAL  PLANE 
As  Defined  by  the  Keplerian  Elements 
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3.3  Models 

Each  of  these  biases  will  be  incorporated  into  the  measurement  model  for  use  in 
this  research.  Biases  often  affect  the  pseudorange  and  phase  range  measurement 
differently,  as  will  be  shown  below. 

Recall  that  this  thesis  is  a  simulation  of  measurements  and  biases.  The  initial, 
imperfect  values  of  the  desired  quantity  will  be  termed  the  "nominal"  value. 
Typically,  corrections  are  applied  to  the  nominal  value  to  reach  the  "true"  value.  In 
Equation  3-17,  the  relationship  between  these  three  quantities  is  shown. 

truth  =  nominal  +  correction  (3-17) 

This  nominal  value  is  the  value  inserted  into  the  measurement  models,  based  on  the 
theory  of  MSOFE  and  its  unique  requirements  for  setting  up  the  problem  at  hand. 
Therefore,  in  all  cases,  the  desired  value  will  be  corrected  as 

desired  (nominal)  =  truth  -  correction  (true  correction)  (3-18) 

where  the  correction  will  be  the  entry  in  the  state  vector,  as  applicable.  As  the  person 
setting  up  the  simulation,  I  shall  supply  the  "truth"  from  which  the  nominal  is  found. 
The  [best  estimate  of  the]  correction  will  be  the  output  of  the  filter. 


3.3.1  Model  for  Pseudorange  Measurements 
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The  final  model  for  the  pseudorange  (in  linear  units)  is  given  by  Equation  3-19, 


T  • 

PR  =  p*  +  c(tR -ts)  +  T*  +-tt  +  S*  sin  <{>  +  C*  cos<|>  +  cAts  +  s  ^  (3-19) 


where 


P*  =  ^[xs -(xR  -8xr)]2  +[ys  ~(yR  -SyR)]2  +[zs  -(zR  -5zr)]2 
=  geometric  distance  between  satellite  at  ts  and  receiver  at  tR 
c(tR  -  ts)  =  speed  of  light  *  clock  differences  due  to  clock  inaccuracies 
where  tR*  =  tR  -5tR 

T*  =  tropospheric  correction  from  Hopfield  model  =  T-6T 

c^/f2  =  speed  of  light  *  ionospheric  correction,  scaled  by  the  square  of  the 

frequency  of  the  carrier  (in  MSOFE,  nominal  =  c(I  -  51)) 

S*  =  arbitrary  coefficient  of  the  sine  term  =  S  -  5S 
C*  =  arbitrary  coefficient  of  the  cosine  term  =  C  -  5C 
4>  ~  location  of  the  satellite  in  the  orbital  plane  based  on  Keplerian  elements 
=  co  +  f 

cAt  =  speed  of  light  *  satellite  clock  correction 
where  Ats  =  Ats  +  Sts 
£  pseudo =  measurement  noise  (uncertainty). 


The  "5"  in  all  cases  represents  the  correction  to  the  nominal 


> 
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3.3.2  Phase  Range  Measurement  Model 


The  final  model  for  the  phase  range  (in  meters)  is  given  by  Equation  3-20, 


cl*  . 

PR  =  p*  +  c(t;  -ts)  +  T*  --V  +  (N xy  +  S*  sin  <J>  +  C*  cos<t>  +  cAts  +  ephasc (3-20) 

f  range 


where 

p*  =  a/[xs  "  (*R  - ' 5xr  )]2  +  [y ! 3  -  (y  R  - ' 5y  R  )]2  +  [zs  -  (zR  -  5zR  jf 
=  geometric  distance  between  satellite  at  ts  and  receiver  at  tR 
c(tR  -ts)  =  speed  of  light  *  clock  differences  due  to  clock  inaccuracies 
where  tR*  =  tR  -6tR 

T*  =  tropospheric  correction  from  Hopfield  model  =  T  -  8T 

-  speed  of  light  *  ionospheric  correction,  scaled  by  the  square  of  the 

frequency  of  the  carrier  (in  MSOFE,  nominal  1  =  c(I  -  51)) 

(NT.)*  =  integer  ambiguity  *  wavelength  of  the  carrier 
S*  =  arbitrary  coefficient  of  the  sine  term  =  S  -  5S 
C*  =  arbitrary  coefficient  of  the  cosine  term  =  C  -  5C 
<}>  =  location  of  the  satellite  in  the  orbital  plane  based  on  Keplerian  elements 
=  co  +  f 

cAt‘  =  speed  of  light  *  satellite  clock  correction 
where  Ats  =  Ats  +  5ts 
ePh«c  =  measurement  noise  (uncertainty) 

range 

The  "5"  in  all  cases  represents  the  correction  to  the  nominal 
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CHAPTER  IV 

RESEARCH  AND  RESULTS 


4.1  Overview 

This  chapter  is  a  summary  of  the  research  leading  to  this  thesis.  Descriptions 
of  the  measurement  and  bias  models  will  be  presented,  as  well  as  the  state  equations 
used  to  define  the  MSOFE  "system"  and  "filter"  models.  The  output  of  MSOFE 
will  show  the  effects  of  different  measurement  and  bias  combinations  on  the  "filter" 
performance  due  to  higher-order  error  (random  bias)  sources  present  in  the  "system 
truth"  model.  Results  will  be  presented  as  plots  of  the  receiver's  position  uncertainty 
as  a  function  of  time 


This  chapter  will  present  the  equations/models  used  by  the  computer  program, 
MSOFE,  for  the  implementation  of  the  Kalman  filter  techniques.  The  problem- 
specific  input  required  by  MSOFE  were1: 


» 


’Recall  that  MSOFE  is  a  "shell"  of  an  optimal  Kalman  filter  processor.  It  consists  of  FORTRAN 
code  divided  up  into  a  main  program  and  84  subroutines.  Only  14  of  these  84  were  required  to  be 
changed  to  meet  the  special  needs  of  the  current  user;  these  user-specific  subroutines  defined  the 
state  equations  and  measurement  models  for  the  system  and  filter.  Other  subroutines  were  written, 
as  needed,  to  augment  MSOFE  and  ensure  all  necessary  information  was  available. 

54’ 
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a.  state  dynamics  (differential)  equations, 

b.  measurement  models, 

c.  partial  derivative  matrices  (used  in  the  linearization  of  the  meas¬ 
urement  models),  and 

d.  the  matrix  to  link  the  true  value  of  the  "filter"  state  vector  to  the 
"system  truth"  state  vector  model  (called  the  "linear  filter/system 
mapping  matrix")  (Carlson  and  Musick,  1990). 

4.2  "System"  versus  "Filter" 

The  distinction  between  the  "system"  and  the  "filter"  must  be  addressed  for  a 
complete  picture  of  the  processes  behind  MSOFE.  As  stated  in  Chapter  1,  within 
MSOFE ,  the  "filter"  is  the  first-order  representation  of  the  real-world,  and  the 
"system"  is  the  second-order  representation.  One  may  also  think  of  the  "filter"  as  an 
under-parameterized  version  of  the  over-parameterized  "system."  This  thesis  set  out 
to  examine  the  effects  on  the  "filter"  performance  due  to  higher-order  error/bias 
sources  in  the  "system  truth"  model.  One  could  think  of  the  "system,"  in  the  context 
of  MSOFE,  as  the  truth  or  real-world  model,  perturbed  by  random  effects  with 
covariance  P.  The  "filter"  is  the  simpler  model  of  the  truth,  containing  a  sufficient 
number  of  states  to  adequately  represent  the  physical  phenomena  of  interest.  The  a 
priori  covariances  within  the  "filter"  will  be  much  greater  than  those  of  the 
"system";  the  purpose  of  the  Kalman  technique  is  to  optimize  the  values  for  the 
desired  quantities  (i.e.,  states  in  the  state  vector)  with  the  help  of  properly  modeled 
measurements. 
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4.3  State  Equation  Models 

The  state  equation  formulation  for  input  into  MSOFE  began  with  Equation  2- 
16,  the  (linearized)  homogeneous  differential  equation  of  the  error  states  with  the 
presence  of  a  noise  input  (recall  from  Section  2.3.2.1  that  the  system  process  noise 
matrix  G  =  I),  seen  in  Equation  4-1. 

6x-F6x  +  w  (4-1) 

Matrix  F  in  Equation  4-1  is  the  homogeneous  state  dynamics  matrix,  6x  is  the 
random  system  error  state  vector,  and  w  is  the  white  noise  vector.  The  GPS  system 
modeled  in  this  thesis  was  modeled  as  a  random  walk  (Gelb,  1992);  therefore. 
Equation  4-1  becomes 


x 


=  w 


(4-2) 


where  F  =  [0],  i.e.,  one  assumes  no  knowledge  about  the  dynamics  of  the  problem. 
A  random  walk  may  be  thought  of  as  a  system  which  makes  fixed-length  steps  in 
arbitrary  directions  (ibid).  Both  the  system  truth  and  the  filter  models  were  modeled 
in  this  fashion.  The  noise  vector  w,  as  Q  =  EfwwT],  was  used  by  MSOFE  in  the 
covariance  propagation  equation  for  the  operational  Kalman  filter,  presented  as 
Equation  2-282.  The  Q  matrix  for  the  "system  truth"  was  the  zero  matrix  because 

2Within  MSOFE,  the  slate  and  covariance  dynamics  equations  are  solved  using  a  single-step 
numerical  integrator  that  uses  a  Runge-Kutta  type  method. 


k  j 
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the  system  was  modeled  as  close  to  ideal  (the  truth)  as  possible  and  did  not  receive 
any  deviations  due  to  insufficient  knowledge  of  the  approximation  process  used  in 
linearization.  The  "filter",  on  the  other  hand,  was  modeled  with  a  process  noise 
matrix.  The  presence  of  this  matrix  indicated  that  the  filter  output  at  time  t+1  did 
not  depend  on  the  output  at  time  t.  In  other  words,  the  presence  of  the  process 
noise  indicated  that  the  time  sequences  of  the  state  vector  and  covariance  contained 
insufficient  information  to  predict  their  respective  values  from  epoch  to  epoch.  Only 
the  states  corresponding  to  the  station  position  and  clock  contained  process  noise, 
implying  that  these  four  states  were  the  unknown  values  to  be  determined  by  the 
filter.  The  "system  truth"  was  modeled  with  extremely  small  a  priori  covariances 
(relative  to  the  filter  model),  as  these  states  are  considered  known 

The  "system"  and  "filter"  error  state  vectors  initially  will  be  identical  (equal 
number  of  states  in  each)  As  biases  are  incorporated  through  the  course  of  the 
research,  the  size  of  the  system  state  error  vector  will  grow,  augmenting  the  "real- 
world"  interpretation  of  the  model.  The  mapping  vector,  mentioned  above,  will  map 
the  like  states  in  the  system  error  vector  to  the  corresponding  true  values  of  the 
filter  states  in  a  one-to-one  manner  (Carlson  and  Musick,  1990),  i.e.,  the  mapping 
matrix  will  be  defined  by  Afs  in  Equation  4-3. 


(4-3) 


58 


Matrix  AfJ  will  be  of  dimensions  (number  of  filter  states)  x  (number  of  system 
states).  This  matrix,  in  this  research,  will  be  defined  as  that  matrix  with  elements 
A(n,n)  =  1,  and  all  other  elements  =  0. 

4.4  Measurement  Models 

As  stated  before,  two  measurement  types  were  used  in  this  research: 
pseudorange  and  phase  range.  The  most  basic  measurement  type  in  GPS  work  is 
the  pseudorange  measurement.  Therefore,  all  iterations  of  the  program  began  with 
the  baseline  using  pseudorange  measurement  data.  To  examine  the  effect  of  adding 
a  more  precise  measurement  type,  the  iterations  were  repeated  using  pseudorange 
and  phase  range  together  through  the  Kalman  filter.  Therefore,  two  measurement 
models  were  necessary.  Table  2  shows  the  measurement  models  for  the 
system/filter  pseudorange  and  phase  range3.  The  symbol  "5"  indicates  an  entry  in 
the  "system"  state  vector  ( true  error),  and  the  "6  hat"  indicates  an  entry  in  the 
"filter"  state  vector  ( estimate  of  the  error).  Measurement  noise,  v,  (uncertainty)  for 
the  pseudorange  measurements  was  set  to  +1  m,  and  measurement  noise  for  the 
phase  range  measurements  was  set  to  +1  mm. 

As  mentioned  previously  in  Chapter  3,  because  of  MSOFE's  unique  approach 
to  Kalman  filtering,  the  "whole  value"  quantity  supplied  by  the  user  is  the  "true" 
value;  the  state  vector  corrects  this  truth  to  become  "nominal"  (linearization  takes 


3  Note  that  the  measurement  biases  were  .added  piecemeal,  and  the  effect  (output)  was  examined  I 

after  each  addition. 
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TABLE  2 

GPS  MEASUREMENT  MODELS 
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ft 


ft 


Pseudorange  Measurement.  "System  Truth"  Model 

ft 

relative  position  from  satellite  (S)  to  receiver 
(R),  R  corrected  for  true  error 

geometric  range  from  S  to  R 

pseudorange  measurement,  time  corrected  for 
true  error 

pseudorange  plus  tropospheric  delay,  ft 

corrected  for  true  error 

pseudorange  plus  ionospheric  delay,  corrected 
for  true  error 

ft 

R4  -  R}+  (A  -  6A)  sin  cp  +  (B  -  bB)  cos  cp  +  c  (AtSAT  -  6  t)  pseudorange  plus 

orbit  biases,  ail  corrected  for  true  error 

R5  ■  R4+  (v„*AvJ  pseudorange  plus  measurement  noise,  scaled  ft 

by  a  normal,  Gaussian  random  number 
generator 


Ir  -Is -(Ir  -SIr) 

P-  I  ll  I 

Ri  =P  +  c((AtR  -btR)-Ats] 
R-,  »  Rj+  Tfscaie  -  6(scale)] 

R3  »  R,+  I[scale  -  6(scale)] 


Pseudorange  Measurement.  "Filter"  Model 

lii  *  rs  -  (r,  -  6r,  +  6r„ )  relative  position  from  S  to  R,  corrected  for 

true  error  and  best  estimate  of  error 

P  =  |  Ir  |  geometric  range  from  S  to  R 


ft 


ft 


ft 
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TABLE  2,  continued 


R,  -p  +  c[(AtR  -6tR  +  6tR)-Ats]  pseudorange  measurement,  time  corrected  for 

true  error  and  best  estimate  of  error 


R2  -  R^  Tfscale] 

pseudorange  plus  tropospheric  delay 

& 

R3  »  R2+  Ifscale] 

pseudorange  plus  ionospheric  delay 

1 

R4  -  R3  +  Asin  cp  +  B  cos  cp  + 

c  AtSAT  pseudorange  plus  orbit  biases 

R5  =  R4  +  vR 

pseudorange  plus  measurement  noise 

» 

"System  Truth”  Model 

&  ”IS"(Ir  -6Ir) 

relative  position  from  S  to  R,  corrected  for 
true  error 

9 

P-  lit  1 

geometric  range  from  S  to  R 

> 

P,  =  p  +  c[(AtR  -StR)-Ats] 

phase  range  measurement,  time  corrected  for 
true  error 

P2  -  P,  +  [NX  -  6(NX)] 

phase  range  plus  integer  ambiguity,  corrected 
for  true  error 

» 

P3  -  P2  +  T[scale  -  6(scale)] 

phase  range  plus  tropospheric  delay, 
corrected  for  true  error 

» 

P4  -  P3  +  I[scaie  -  6(scale)] 

phase  range  plus  ionospheric  delay,  corrected 
for  true  error 

P5  -  P4  +  (A  -  6A)  sin  cp  +  (B  - 

6B)  cos  cp  +  c  (AtSAT  -  5  t)  phase  range  plus  orbit 
biases,  all  corrected  for  true  error 

» 

» 


TABLE  2,  continued 


P6  -  P5  +  (v,  *Av,)  phase  range  plus  measurement  noise,  scaled 

by  a  normal,  Gaussian  random  number 
generator 


Phase  Ranee  Measurement.  "Filter”  Model 


Is,  -Is  -(I,  -6r„  +6r() 


relative  position  from  S  to  R,  corrected  for 
true  error  and  best  estimate  of  error 


P  =  Is 


geometric  range  from  S  to  R 


P.  -  p  +  c[(AtR  -  6tR  +  &tR) -  Ats  ]  phase  range  measurement,  time  corrected  for 

true  error  and  best  estimate  of  error 

P,  -  P,+  [NX  -  6(NX)  +  6(NX)J  phase  range  plus  integer  ambiguity,  corrected 

for  true  error  and  best  estimate  of  error 


P3  -  P,  +  T[scale] 


phase  range  plus  tropospheric  delay 


P4  .  p.  +  Ifscale] 


phase  range  plus  ionospheric  delay 


P5  -  P4  +  A  sin  cp  +  B  cos  cp  +  c  At  phase  range  plus  orbit  biases 


P„-P5+  v. 


phase  range  plus  measurement  noise 
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place  around  the  nominal  value).  Recall  Equation  3-18  which  shows  the  relationship 
between  the  "truth",  the  "nominal",  and  the  "correction".  Using  this  logic,  the 
"system"  measurement  model  in  MSOFE  will  correct  the  desired  quantities  as  in 
Equation  4-4,  and  the  "filter"  measurement  model  as  in  Equation  4-5. 

SYSTEM:  nominal  =  truth  (whole  value)  -  true  error  (4-4) 

FILTER:  nominal  =  truth  -  true  error  +  best  estimate  of  the  error  (4-5) 

Data  used  in  the  measurement  models  were  simulated.  All  efforts  were  made 
to  recreate  as  realistic  a  measurement  value  as  possible.  Twenty-eight  possible 
satellite  numbers  were  examined  for  each  measurement  type,  for  a  total  of  56 
possible  measurements.  Obviously,  all  28  satellite  numbers  did  not  contribute  to  the 
update,  as  information  on  every  satellite  was  not  included  in  the  broadcast 
ephemeris  Not  knowing  in  advance  which  satellites  were  listed  in  the  ephemeris, 
provisions  were  made  for  all  possible  satellite  combinations.  Briefly  summarized, 
here  is  a  list  of  steps  taken  to  generate  the  necessary  information  with  which  the 
range  data  were  created  within  the  software: 

a  Broadcast  ephemeris  file  was  read  [by  the  computer  program]  to 
calculate  the  nominal  satellite  position  coordinates  (per  epoch), 
ionospheric  correction  coefficients,  and  the  location  of  the 
satellites  within  the  orbital  plane  (argument  of  latitude,  4>). 
b  Receiver  data  file  for  the  corresponding  broadcast  ephemeris 
was  read  to  return  the  nominal  station  coordinates. 
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c.  Elevation  of  each  satellite  was  calculated  (per  epoch). 

d.  Tropospheric  delay  of  each  satellite  was  calculated  (per  epoch), 
based  on  the  elevation. 

e.  Ionospheric  delay  of  each  satellite  was  calculated  (per  epoch), 
based  on  the  elevation. 

f.  Orbital  errors  for  each  satellite  were  calculated  (per  epoch), 
based  on  <j>. 


~ — 
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ft 


ft 
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ft 
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A  point  must  be  made  concerning  the  epochs  used  in  the  research.  Time  tags 
of  each  epoch  were  incremented  by  MSOFE  in  order  to  update  the  state  vectors  and 
covariances  Updates  were  made  every  "second"  (from  0  to  100  seconds),  meaning 
that  each  "second,"  measurements  were  taken,  at  that  time,  an  update  was  made. 
Time  tags  used  to  generate  the  measurements  were  updated  each  time  a 
measurement  was  needed.  In  order  to  see  a  distinct  movement  of  the  satellites  from 
the  ephemeris  file  information,  the  time  sent  to  the  computer  program's  subroutines 
to  calculate  the  measurements  was  incremented  by  30  "seconds."  In  other  words, 
the  time  used  by  MSOFE  and  that  used  by  the  measurement-generation  routine  were 
different,  but  they  were  each  transparent  to  one  another.  The  important  aspect  of 
the  time  definitions  was  to  ensure  that  MSOFE  and  the  measurement-generation 
routine  conducted  updates  at  the  same  instant— their  relative  times  were,  as  such. 


ft 


ft 


unimportant. 


In  all  cases,  the  station  coordinates  and  the  station  clock  correction  were 
included  in  the  filter  and  system  state  vectors.  The  a  priori  value  for  the  standard 
deviations  (in  meters)  of  these  quantities  may  be  seen  in  Equation  4-6  and  4-7. 

[  aSx  a6y  a6z  CTSt  U  = 

[±lxlO“*±Ixl(r*  ±1x1  (T4  ±1x1  (T4]7  (4-6) 

[  °Sx  CTSv  CTSz  °5t  ]f,|ter  = 

[±lxlO_s±lxlO-2  ±lxlO“2  ±lxl(T2 ]T  (4-7) 

For  iterations  containing  phase  range  data,  the  integer  ambiguity  situation  had 
to  be  addressed.  The  integer  was  considered  an  arbitrary  number,  so  any  value  for 
an  initial  integer  would  have  sufficed;  an  arbitrary  value  of  lxlO6  was  chosen  for 
each  satellite.  The  a  priori  standard  deviation  (in  meters)  of  the  integer  ambiguity 
was  chosen  as  ±1x10^  for  the  system  and  ±1x1 02  for  the  filter.  This  bias  was 
present  in  both  the  "system"  and  "filter"  state  vectors,  but  only  when  phase 
measurements  were  processed. 

The  model  for  the  tropospheric  delay  was  taken  from  the  Hopfield  model 
(Goad  and  Goodman,  1974).  The  model  creates  the  inverse  relationship  between 
the  satellite's  elevation  and  the  correction  due  to  the  troposphere— the  higher  the 
satellite  with  respect  to  the  station,  the  less  the  delay  due  to  the  troposphere.  This 
inverse  relationship  is  logical  when  one  thinks  of  how  the  troposphere  spans  away 
from  the  earth's  surface.  If  the  satellite  were  high  over  the  station,  the  GPS  signal 
would  travel  less  through  the  propagation  medium  to  reach  the  receiver  than  if  the 
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satellite  were  located  on  the  horizon.  The  scale  for  the  delay  was  :hosen  as  1.0. 
The  entry  in  the  state  vector  contained  the  uncertainty  of  the  scale,  as  shown  in 
Table  2.  The  a  priori  value  for  the  standard  deviation  of  the  scale  in  the  system 
model  was  chosen  as  +2  percent4  (ibid).  This  bias  was  not  included  in  the  "filter" 
model. 

The  model  for  the  ionosphere  was  taken  from  the  ICD  for  the  standard  GPS 
receiver  (Leick,  1990).  This  model,  too,  shows  the  inverse  relationship  between  the 
elevation  of  the  satellite  with  respect  to  the  receiver  and  the  amount  of  ionospheric 
delay.  The  scale  for  the  delay  was  chosen  as  1.0.  The  entry  in  the  state  vector 
contained  the  uncertainty  of  the  scale,  as  shown  in  Table  2.  The  a  priori  value  for 
the  standard  deviation  of  the  scale  in  the  system  model  was  chosen  as  +50  percent 
(ibid)  This  bias  was  not  included  in  the  filter  model 
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The  model  for  the  orbit  biases  was  taken  from  current  theory  of  the 
perpetuation  of  errors  in  the  orbit.  This  model  constitutes  a  "once-per-revolution" 
effect  of  the  perturbations  of  the  Keplerian  elements,  which  define  the  orbit  of  the 
satellite.  The  periodic  oscillations  of  the  satellite  as  it  travels  along  its  orbit  and  the 
deviation  of  the  satellite  clock  are  modeled  in  Equation  4-8.  This  model  represents 
the  change  in  the  satellite-receiver  range  due  to  orbital  biases.  The  argument  of 
latitude  ((j>)  places  each  satellite  in  its  respective  orbit.  The  clock  term  represents  the 
deviation  of  the  satellite  clock  from  predicted  models  in  the  broadcast  ephemeris 


4Found  by  rms  error  of  the  Hopfield  model  divided  by  the  average  correction. 


66 


Ar  =  A  sin  cj>  +  B  cos<t>  +  At"  (4-8) 

information.  Coefficients  A  and  B  represent  the  strength  of  each  oscillation.  In 
MSOFE,  the  coefficients  were  set  to  zero;  the  a  priori  value  for  the  standard 
deviation  of  covariance  of  the  coefficients  was  10  m  (initially,  but  other  iterations 
were  made  with  a  deviation  of  + 1  m,  for  comparison)  The  standard  deviation  of 
the  satellite  clock  bias  was  initially  set  to  + 1  ns.  These  biases  were  only  present  in 
the  "system"  state  vector. 

As  stated  above,  the  computer  program  allowed  for  as  many  as  28  different 
satellite  numbers  Therefore,  28  different  measurements  were  possible  when 
processing  pseudorange  data,  and  28  were  possible  when  processing  phase  range 
data.  For  the  satellite-specific  biases  (integer  ambiguity  and  the  orbital  biases),  the 
state  vector  had  to  account  for  one  state  per  possible  satellite.  Table  3  shows  the 
size  of  the  "filter"  and  "system"  state  vectors  as  the  number  of  measurement  biases 
changed  The  integer  ambiguity  bias  was  not  present  in  the  state  vectors  for 
iterations  using  pseudorange  data  only. 

The  last  condition  investigated  was  that  of  cycle  slips.  As  mentioned  before,  a 
cycle  slip  occurs  wuen  the  signal  from  satellite  to  receiver  is— for  whatever  reason- 
disrupted.  This  disruption  often  stems  from  an  actual,  physical  blocking  of  the 
signal  from,  for  example,  a  tree  or  building,  or  perhaps  the  GPS  receiver  on  an 
aircraft  is  blocked  from  the  signal's  path  by  the  plane's  own  maneuvering.  In  the 
filter,  the  cycle  slip  was  simulated  by  setting  up  a  flag  to  look  for  Satellite  #3  during 
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epochs  3,  20,  and  70.  When  these  times  occurred,  a  large  value  for  Q,  the  process 
measurement  noise,  was  inserted  in  the  program  associated  with  Satellite  #3.  The 
large  noise  value  told  the  program  that  this  value  was  now  unknown  (at  those  state 
epochs),  and  the  effect  of  the  slip  was  analyzed  on  output  of  the  filter. 


ft 


ft 


* 


* 


» 


TABLE  3 

STATE  VECTOR  ENTRIES  PER  INPUT 


BIAS 

SYSTEM  VECTOR 

FILTER  VECTOR 

x  coordinate 

1 

1 

y  coordinate 

1 

1 

z  coordinate 

1 

1 

receiver  clock 

l 

1 

28 

28 

| 

1 

not  present 

ionosphere 

1 

not  present 

orbit  (sine  term) 

28 

not  present 

orbit  (cosine  term) 

28 

not  present 

satellite  clock 

28 

no!  present 
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4.5  Results 

This  section  will  present  the  results  of  the  iterations  through  the  Kalman  filter 
program  as  measurement  type  and  biases  were  varied.  The  presentation  of  the 
output  will  be  categorized  by  bias  type.  Then,  comparisons  will  be  made  between 
the  pseudorange  and  pseudo/phase  range  combinations.  The  different  system  state 
vector  input  combinations  for  each  output  graph  may  be  seen  in  Tables  4,  and  55. 
The  reader  is  referred  to  the  end  of  this  chapter  for  all  output  figures. 

4.5.1  Geometry  and  Dynamics 

An  important  aspect  of  GPS  theory  is  the  relationship  of  the  satellites  to  the 
respective  receiver.  The  geometry  of  one  with  respect  to  the  other  plays  a  major 
role  in  the  user's  ability  to  recover  the  receiver's  position  Recall  that  the  broadcast 
ephemeris  was  incremented  in  time  to  simulate  the  movement  of  the  satellites.  Two- 
dimensional  plots  of  the  satellites'  location  relative  to  the  receiver  may  be  seen  in 
Figures  8,  9,  10,  and  1 16.  The  figures  show  the  satellites'  relative  locations  at  four 
instances  in  the  time  interval  of  study.  Without  the  presence  of  any  other  anomalies, 
the  relative  spacing  of  the  satellites  at  times  1  and  60  seconds  (Figures  8  and  9, 
respectively)  would  indicate  a  better  possibility  for  recovering  one's  position  with 
greater  certainty.  Further,  the  presence  of  a  fifth  satellite  at  t=60  sec  further 

5Not  all  combinations  of  biases  and  measurement  types  are  presented  here.  Often,  changing  a 
variable  did  not  appreciably  change  the  output--these  plots  were  not  included  because  they  did  not 
offer  any  significant  findings 

6Station  and  satellite  coordinates  were  presented  in  the  topocentric  system. 
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orbit  coefficient  =  I  m 

error  on  one  satellite  only 

elev-dependent  (tropo) 

elev-dependent  (ionos) 

equal  msmt  uncertainties 

cycle  slip 

SYSTEM  SET-UP 

MEASUREMENT 

TYPE 

&  DYNAMICS 

STATION- 

DEPENDENT  BIASES 

OBSERVATION- 
DEPENDENT  BIASES 

SATELLITE- 
DEPENDENT  BIASES 

MISCELLANEOUS 

increases  the  station/satellite  polyhedron's  "volume"  (mentioned  in  Chapter  3), 
possibly  leading  to  an  even  better  position  recovery.  Figures  1 0  and  1 1  show  the 
geometry  begin  to  degrade  as  the  satellites  move  topocentrically  away  from  the 
receiver.  One  would  expect,  in  an  ideal  situation,  the  position  uncertainty  to 
increase  as  the  geometry  degrades. 

In  addition  to  the  geometry  of  the  satellites  and  the  receiver,  the  elevation  of 
the  satellite  with  respect  to  the  receiver  is  also  an  important  issue.  The  reader  is 
referred  to  Figures  12  and  13,  plots  of  the  satellites'  elevation  and  elevation  versus 
azimuth,  respectively,  over  the  study  interval.  This  plot  was  generated  using  the 
information  from  the  broadcast  ephemeris7.  As  stated  above,  the  amount  the  GPS 
signal  is  affected  by  the  propagation  media  is  inversely  proportional  to  the  sine  of 
the  satellite’s  elevation.  One  may  expect  the  position  uncertainty  to  increase  as  the 
elevation  angles  decrease  due  to  the  larger  bias  affecting  the  measurement. 

The  data  reduction  through  MSOFE  began  with  the  non-practical  case  of  a 
non-moving  receiver  with  non-moving  satellites.  In  this  case,  station  and  satellite 
approximate  coordinates  were  hard-programmed  into  MSOFE  (were  not  taken  from 
the  ephemeris,  as  on  subsequent  runs).  The  purpose  was  to  examine  the  uncertainty 
of  the  station's  position  when  pseudorange  measurements  only  were  processed,  with 
no  external  influences.  Next,  the  consequence  of  adding  phase  range  measurements 
was  examined.  Figure  14  shows  the  results.  The  plot  of  pseudorange  only  shows  a 


7The  satellite  numbers  shown  are  those  with  data  present  in  the  ephemeris  whose  elevation  was 
greater  than  a  pre-defmed  cut-ofT  value  of  0  degrees. 
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steady  uncertainty,  logical  because  the  station  and  satellite  were  stationary.  The 
addition  of  the  phase  ranges  immediately  improves  the  overall  uncertainty  of  the 
position;  the  phase  bias  for  each  satellite  contributes  a  standard  deviation  at  the  rate 
of  V 'j-  where  n  is  the  number  of  measurements  made. 


In  Figure  15,  the  station  and  satellite  coordinates  were  taken  from  the 
broadcast  ephemeris,  varied  with  time.  Due  to  the  degradation  of  the  geometry,  the 
uncertainty  of  position  using  the  pseudoranges-only  model  increases  from  t=0  to 
about  t=83.  The  slight  "glitch"  at  t=58  is  due  to  the  addition  of  the  fifth  satellite. 
The  uncertainty  of  position  decreases  for  (1-2  seconds)  a  small  amount  due  to  the 
added  redundancy8,  but  the  low  elevation  of  this  fifth  satellite  and  the  further 
degradation  of  the  geometry  is  counter-productive  in  keeping  the  uncertainty  to 
within  acceptable  levels.  The  incorporation  of  the  phase  range  clearly  enhances  the 
system,  and  the  position  uncertainty  at  the  end  of  the  study  time  approaches  the 
millimeter  level.  The  continued  rising  of  satellites  #12  and  16  (see  Figure  12)  past 
t=83  helps  to  further  reduce  the  uncertainties. 

4.5.2  Propagation  Media 


* 


i 
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The  next  stage  of  study  was  to  add  the  effects  of  the  propagation  media  Four 
iterations  were  made  for  each  type  of  measurement— with  and  without  tropospheric 


8Recall  that  the  quantities  being  estimated  are  the  station  coordinates  and  the  station  clock 
correction.  A  minimum  of  four  satellites  are  required  to  solve  for  the  four  unknowns.  Any 
satellites  beyond  the  required  four  create  a  redundant  situation  and  should  decrease  uncertainties, 
providing  all  other  variables  (elevation,  propagation  media  delay,  etc.)  are  also  acceptable. 
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biases,  and  with  and  without  ionospheric  biases.  Table  6  shows  the  four 
combinations  and  their  other  names  that  will  be  used  interchangeably  throughout 
this  thesis. 
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TABLE  6 

PROPAGATION  MEDIA  COMBINATIONS 


COMBINATIONS 

ALTERNATE  NAME 

ERROR-FREE  TROPOSPHERE 
ERROR-FREE  IONOSPHERE 

IDEAL  CASE 

(IDEAL/  DUAL-FREQUENCY ) 

ERROR  ON  TROPOSPHERE 
ERROR-FREE  IONOSPHERE 

DUAL-FREQUENCY  CASE 

ERROR-FREE  TROPOSPHERE 
ERROR  ON  IONOSPHERE 

<  NONE  > 

ERROR  ON  TROPOSPHERE 
ERROR  ON  IONOSPITERE 

REAL-WORLD  CASE 
(SINGLE-FREQUENCY  CASE) 

The  "ideal"  case  contains  no  propagation  media  errors.  The  dual-frequency  case  is 

so  named  because  with  a  dual-frequency  GPS  receiver,  it  is  possible  to  eliminate  the  > 

linear  affects  of  the  ionosphere  using  post-processing  measurement  models  (due  to 

the  dispersive  characteristics  of  the  ionosphere).  The  third  combination  is  not 

practical  in  the  real  world— the  effect  due  to  troposphere  cannot  be  eliminated  while  I 


> 


74 


the  ionospheric  effect  remains.  This  combination  was  produced  for  comparison 
purposes  only  and  will  not  be  discussed  in  detail9.  The  "real-world"  case  is  so 
named  because  it  contains  errors  on  both  media.  It  is  also  known  as  the  single¬ 
frequency  case  because  both  biases  will  be  present  with  a  single-frequency  receiver 
(neither  can  be  eliminated  without  further  information). 

Figure  16  shows  the  effect  of  the  propagation  media  on  the  model  containing 
only  pseudorange  measurements  The  ideal  case  (no  errors  on  either)  is  the  same 
plot  as  seen  in  Figure  15.  The  effect  of  adding  the  troposphere,  known  to  2%,  to 
get  the  dual-frequency  model  only  increases  the  uncertainty  by  less  than  0.1  m, 
shown  in  Figure  1710.  However,  adding  the  ionospheric  bias  with  a  standard 
deviation  of  50%  increases  the  uncertainty  by  approximately  2.5  m.  In  the  single¬ 
frequency  plot,  the  "glitch"  at  t=58  sec  is  more  pronounced  than  that  in  the  plot  void 
of  these  biases,  more  clearly  seen  in  Figure  18.  Recall  at  t=58  that  the  elevation  of 
the  fifth  satellite  rises  to  above  0  degrees  and  enters  into  the  filtering  process.  The 
low  elevation  contributes  a  large  amount  of  delay  to  the  measurement  (about  84  m 
tropospheric  correction  and  21  m  ionospheric  correction).  The  redundancy  of  the 
additional  satellite  is  not  enough  to  counteract  the  large  contribution  of  the 
propagation  media,  so  the  uncertainty  rises  by  approximately  2  m  Contrast  this 
jump  in  uncertainty  to  the  slight  decrease  of  uncertainty  in  the  ideal  case,  where 

9This  combination  typically  showed  little  difference  from  the  ideal  case  due  to  the  small  (  +  2%) 
uncertainty  on  the  tropospheric  delay. 

l0NOTE:  The  graphs  with  and  without  the  tropospheric  bias  (with  the  same  ionospheric  biases) 
were  essentially  the  same.  The  plots  were  offset  slightly  on  the  plot  so  that  the  four  plots  could  be 
distinguished  from  one  another. 
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elevation  does  not  play  a  role  in  determining  the  uncertainty.  Subsequent  to  t=58, 
the  degrading  geometry  somewhat  blends  with  the  effect  of  elevation,  i.e ,  satellites 
#3,  17,  and  26  are  sufficiently  high  above  the  receiver  so  to  not  degrade  the 
uncertainty,  despite  the  geometry. 

The  addition  of  phase  ranges  to  the  pseudorange  measurements  may  be  seen  in 
Figure  19  for  all  four  combinations  of  propagation  biases.  The  ideal  case  is  the 
same  plot  as  seen  in  Figure  15.  Despite  the  mm-level  measurement  noise  of  the 
phase  range  measurements,  the  addition  of  the  tropospheric  bias  creates  a 
pronounced  increase  in  the  position  uncertainty.  From  t=0  to  t=58  sec,  the  ideal 
case  and  the  dual-frequency  case  are  identical  at  the  millimeter  level  (see  Figure  20 
for  a  clearer  look  at  the  two  graphs)  With  the  addition  of  the  fifth  satellite,  the 
redundancy  cannot  make  up  for  the  low  elevation  and  the  associated  large 
uncertainties  imposed  on  the  measurements.  Figure  21  shows  the  ideal  case  versus 
the  real-world  case  The  addition  of  tropospheric  and  ionospheric  biases  to  the 
combined  pseudo/phase  range  model  greatly  increases  the  uncertainty  of  the 
receiver's  position  The  redundancy  of  one  satellite  cannot  correct  for  the  large 
contribution  the  satellites'  elevations  pay  toward  the  uncertainty.  The  three  setting 
satellites  and  the  low  (redundant)  satellite  contribute  too  much  uncertainty  to  the 
problem. 

Figure  22  provides  a  look  at  the  comparison  between  the  two  measurement 
types.  Between  t=0  and  t=58  sec,  the  pseudo/phase  range  combination  provides 
better  results  than  the  pseudorange-only  model  Affer  t=58  sec,  the  pseudorange- 
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only  case  is  comparatively  unaffected  by  the  poor  geometry  and  low  elevation  angle 
of  the  redundant  satellite.  However,  the  pseudo/phase  range  model  shows  a  greater 
reaction  to  the  elevation-dependent  biases.  Recall  from  Chapter  3  that  the 
ionosphere  affects  the  phase  and  pseudorange  values  in  opposite  directions.  The 
combination  of  the  two  measurement  types  actually  magnifies  the  effect  of  the  bias 
on  the  recovered  position  uncertainty.  As  the  fifth  satellite  rises  to  a  more 
acceptable  level  (around  10  degrees  elevation),  and  the  maintenance  of  three  other 
satellites  at  acceptable  elevations,  the  pseudo/phase  range  and  the  pseudorange-only 
cases  approach  the  same  level  of  uncertainty. 

4.5.3  Elevation-Dependency 

Due  to  the  results  of  the  propagation  media  models  in  the  previous  section, 
another  model  was  developed— the  elevation-dependent  model.  Two  models  were 
developed:  one  based  on  the  tropospheric  correction  values  and  one  based  on  the 
ionospheric  correction  values.  In  these  models,  the  measurement  noise  was 
weighted  to  give  less  emphasis  to  those  satellites  with  low  elevation  angles  and  more 
emphasis  to  those  with  higher  elevation  angles.  The  new  uncertainty  may  be  seen  in 
Equations  4-9  and  4-10  for  the  tropospheric-dependent  and  ionospheric-dependent 
models,  respectively, 

=  V(°™.)2+(0«P.)2  =  ):+( 0.02T)2  (4-9) 

=  V(^,)2+(  =  V(On-n,)2+(0.5I)2  (4-10) 
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where  T  and  I  are  the  calculated  values  for  the  tropospheric  and  ionospheric 
corrections,  respectively,  and  the  omsml  is  the  standard  deviation  of  the  measurement 

(phase  range  or  pseudorange).  Using  these  models,  one  can  see  that  the  combined 
deviation  will  be  greater  than  that  of  just  the  measurement  deviation.  The  standard 
deviation  of  the  measurement  noise  will  be  used  to  calculate  the  filter  updates  per 
the  Kalman  update  equations  and  will  have  an  overall  effect  of  R“'  (inverse  of  the 
measurement  noise  matrix,  see  Chapter  2)  on  the  update.  Therefore,  satellites  with 
large  tropospheric/ionospheric  corrections  will  contribute  less  to  the  process. 

Because  the  problem  requires  four  satellites  to  determine  uniquely  the  four 
unknowns,  no  change  in  the  output  of  the  filter  should  be  expected  when  only  four 
satellites  are  present.  Therefore,  the  weightings  will  only  come  into  play  with  the 
addition  of  the  redundant  satellite(s).  One  would  expect,  therefore,  to  see  little 
change  in  the  uncertainty  curves  (with  and  without  the  model)  up  to  t=58  sec,  and 
the  addition  of  Satellite  #12  beyond  t=58  should  improve  the  uncertainty  results. 

4.5.3. 1  Dependency  Based  on  Troposphere 

Figure  23  shows  the  pseudorange-only  uncertainty  curves  with  the 
measurement  noise  dependent  on  the  tropospheric  corrections.  As  expected,  the 
uncertainty  curves  for  the  real-world  cases  with  and  without  the  elevation 
dependency  do  not  differ  when  only  four  satellites  are  visible  above  the  horizon. 
However,  with  the  addition  of  the  fifth  satellite,  the  position  uncertainty  does 


decrease,  albeit  only  slightly.  The  small  uncertainty  of  the  troposphere  correction 
(only  ±2%)  does  not  play  an  extensive  role  in  the  weighting  of  the  measurement 
noise,  and  not  much  improvement  is  noted. 

Figure  24  shows  the  pseudo/phase  range  uncertainty  curves  with  the 
measurement  noise  dependent  on  the  tropospheric  corrections.  As  expected,  there 
is  no  change  in  the  uncertainty  curves  with  and  without  the  elevation-dependent 
model  through  t=58  sec.  Dramatic  results  occur,  however,  after  t=58.  The  large 
spike  in  the  uncertainty  after  the  addition  of  the  fifth  satellite  does  not  occur  in  the 
elevation-dependent  model.  Instead,  the  uncertainty  curve  is  smooth;  the  weighting 
of  the  lowest  (in  elevation)  satellite's  measurement  plays  down  the  importance  of 
this  measurement.  The  dip  in  the  elevation-dependent  curves  at  approximately  t=86 
sec  occurs  near  the  point  where  Satellites  #12  and  #16  reach  equal  elevations.  The 
weighting  tends  to  create  a  non-optimal  geometry  by  de-emphasizing  the 
measurements  from  the  satellites  with  the  lowest  elevation  Past  t=86  sec,  the 
measurement  from  Satellite  #12  is  de-emphasized.  If  one  recalls  Figure  10,  the 
geometry  of  the  station  and  satellites,  one  can  see  the  effect  on  the  geometry  of 
(essentially)  removing  Satellite  #12.  The  "polyhedron"  is  greatly  reduced  in  size, 
thereby  increasing  the  position  uncertainty.  Past  t«95  sec,  the  measurement  for 
Satellite  #16  is  de-emphasized;  one  can  see  from  Figure  11  the  effect  on  the 
geometry  of  eliminating  Satellite  #16  at  this  time.  Here,  the  implementation  of  the 
elevation-dependent  measurement  noise  model  does  not  produce  optimum 
uncertainty  results.  The  three  setting  satellites  and  the  degrading  geometry  cannot 
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be  compensated  for  in  the  elevation-dependent  model,  and  the  elevation-independent 
model  yields  more  realistic  results. 

4.5.3.2  Dependency  Based  on  Ionosphere 

Figure  25  shows  the  elevation-dependent  measurement  noise  model,  based  on 
the  ionospheric  correction,  in  the  pseudorange-only  measurement.  Just  as  for  the 
troposphere-dependent  model  in  Figure  23,  the  elevation-dependent  and 
independent  models  are  identical  prior  to  t=58  sec.  Subsequent  to  t=58,  in  the 
presence  of  the  redundant  satellite,  the  uncertainty  increases.  Clearly,  there  is  a 
trade-off  between  weighting  of  the  measurements  and  the  uncertainty  of  the 
position.  By  giving  less  weight  to  the  lowest  (elevation)  satellite,  the  geometry  is 
degraded,  and  the  uncertainty  increases. 

In  Figure  26,  the  elevation-dependency  is  based  on  the  ionospheric  correction 
model,  and  the  measurement  model  is  the  pseudo/phase  range  combination.  Here, 
the  presence  of  the  phase  range  measurements  slightly  increases  the  uncertainty  (less 
than  1  m)  for  epochs  prior  to  t=58  sec.  However,  after  t=58  sec,  the  elevation- 
dependent  model  yields  improved  results.  There  is  still  some  trade-off"  in  the 
geometry  at  epochs  greater  than  97  sec,  but  overall,  the  uncertainty  is  improved 
with  the  elevation-dependent  model.  The  addition  of  phase  ranges  adds  to  the 
improvement  over  the  case  seen  in  Figure  25. 
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4.5.4  Orbital  Biases 

The  orbital  errors  were  added  in  two  different  ways.  First,  the  trigonometric 
coefficients  were  given  an  a  priori  standard  deviation  of  +10  m.  Second,  the 
trigoncmetiic  coefficients  were  given  a  standard  deviation  of  + 1  m.  In  both  cases, 
the  satellite  clock  was  given  an  a  priori  standard  deviation  of  + 1  ns.  Each  way  will 
be  addressed. 
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Figure  27  shows  the  pseudorange-only  case  for  the  o  =  +10  m  case.  The 
measurement  models  with  and  without  the  orbital  errors  are  contrasted  The  figure 
clearly  shows  that  the  large  uncertainty  of  the  orbit  contributes  to  the  increased 
uncertainty  of  the  position  of  the  station.  The  increase  in  uncertainty  is 
approximately  ten-fold  at  its  peak,  at  approximately  t=83  sec.  In  Figure  28,  the 
standard  deviation  of  the  orbital  bias  coefficient  is  set  to  +  1  m  Obviously,  knowing 
the  coefficients  ten  times  better  serves  to  reduce  the  uncertainty  in  the  position  by 
approximately  the  same  amount  at  its  peak  (around  t=83  sec). 
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Figure  29  shows  the  real-world  model  using  pseudoranges  only  and  an 
elevation-dependent  measurement  noise  model,  based  on  the  ionosphere.  Curves 
with  both  values  for  and  without  the  orbital  bias  are  shown.  The  addition  of  the 
elevation-dependency  mode!  in  concert  with  the  orbital  bias  problem  did  not  change 
the  plots  from  those  seen  in  Figures  27  and  28.  The  presence  of  the  orbit  error  was 
significant  enough  to  override  the  weighting  of  the  lower  satellite.  No  significant 
value  was  gained  by  using  the  elevation-dependency  model  in  this  case. 
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In  Figure  30,  the  addition  of  phase  ranges  to  the  "filter"  observation  model 
vastly  improves  the  uncertainty  over  the  pseudorange-only  case  of  Figure  27.  In 
Figure  30,  the  orbit  bias  using  a  coefficient  a  =  +10  m,  through  t=58  sec,  is 
identical  to  the  uncertainty  curve  of  Figure  27  (pseudorange  only).  With  the 
addition  of  the  redundant  satellite,  the  uncertainty  drops  markedly  and  follows  the 
same  shape  as  the  model  without  the  orbital  errors.  In  contrast  to  the  curve  void  of 
orbital  biases,  the  uncertainty  curve  with  the  orbital  bias  differs  by  only  as  much  as  9 
m  at  its  largest  point,  a  stark  improvement  from  the  comparison  shown  in  Figure  27. 
Although  the  orbital  bias  adds  uncertainty  to  the  position  of  the  station,  the  addition 
of  phase  ranges  improves  the  results  overall. 

In  Figure  3 1,  the  curves  found  using  o  =  + 1  m  are  all  virtually  identical  prior 
to  t=58  sec.  The  addition  of  the  orbit  errors  subsequent  to  that  time  yields  results 
which  parallel  those  of  the  no-orbit-bias  curve.  The  presence  of  the  redundant 
satellite  does  not  counteract  the  presence  of  the  tropospheric,  ionospheric,  and 
orbital  biases. 

Figure  32  looks  at  the  situations  already  shown  in  Figures  30  and  31,  except 
with  the  presence  of  the  elevation-dependent  measurement  noise  model  based  on  the 
ionosphere.  The  combination  of  the  uncertainty  in  the  ionosphere  together  with  the 
ionospheric,  tropospheric,  and  orbital  biases  creates  a  situation  which  produces  an 
unacceptable  uncertainty  curve  when  compared  to  other  schemes. 
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Figure  33  investigates  the  contribution  of  one  satellite,  here  Satellite  #3, 
infected  with  orbital  biases,  where  the  orbit  coefficients'  a  priori  ct  =  +10  m  for 
Satellite  #3,  and  a  =  +lxlO~20  m  for  all  other  satellites.  The  overall  uncertainty 
curve  found  from  a  biased  Satellite  #3  is  seen  to  contribute  approximately  10  m  to 
the  case  without  orbital  biases. 

4.5.5  Equal  Measurement  Noise 
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The  next  situation  examined  was  that  of  equal  measurement  noise.  Both  the 
pseudorange  and  phase  range  measurements  were  assigned  a  measurement  noise  of 
+ 1  m  in  the  filter  model.  Figure  34  may  be  compared  to  Figure  21.  In  Figure  34, 
the  presence  of  the  equal  measurement  noise  does  not  improve  the  uncertainty  from 
the  beginning  for  the  ideal  case.  For  the  single-frequency  model,  the  curve  in  Figure 
21  and  Figure  32  are  identical  up  until  t=58  sec.  Although  the  "spike"  from  Figure 
21  does  not  appear  after  that  time,  the  equal-noise  case  does  not  lead  to  the 
improved  uncertainty  in  the  face  of  the  degraded  geometry  as  time  approaches 
t=l  00  sec. 

In  Figure  35,  the  elevation-dependent  (based  on  the  ionosphere)  model  is 
examined.  There  is  definitely  a  trade-off  past  t=58  sec  in  the  uncertainty  curves  for 
the  elevation-dependent  and  -independent  models.  As  the  geometry  degrades,  the 
uncertainty  improves  with  the  former  model;  the  latter  model  sees  the  better  results 
as  the  redundant  satellite  rises  in  the  sky,  but  before  the  geometry  degrades.  As  in 
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Figure  34,  the  equal  measurement  noises  do  not  serve  to  improve  the  overall  output 
of  the  filter,  hence  the  uncertainty  of  the  station's  position. 

In  Figures  36  and  37,  the  orbital  errors  are  introduced  (c  =  +10  and  +1  m, 
respectively).  The  uncertainty  curve  resembles  those  in  Figures  27  and  28, 
respectively,  for  the  pseudorange-only  case.  The  uncertainty  in  Figure  36  at  its  peak 
is  slightly  better  than  in  Figure  27,  possibly  due  to  the  increased  number  of 
measurements  in  the  pseudo/phase  range  situation  and  the  differences  in  the 
measurement  models  for  each  case.  The  same  case  could  be  made  for  Figures  37 
and  28. 

Figure  38  shows  the  case  in  Figure  36  with  the  elevation-dependent  (based  on 
ionosphere)  model.  This  curve  resembles  that  of  Figure  32,  but  this  curve  peaks  at  a 
smaller  uncertainty  value,  for  the  reasons  mentioned  above. 

4.5.6  Cycle  Slip 

One  final  iteration  was  made  for  this  problem.  The  issue  of  cycle  slips,  or 
losses  of  lock,  was  analyzed.  Figure  39  shows  the  results  of  placing  a  simulated 
cycle  slip  at  t  =  3,  20,  and  70  seconds.  When  only  four  satellites  are  visible,  the 
cycle  slip  has  no  affect  on  the  uncertainty  curve.  By  introducing  a  cycle  slip,  the 
program  is  being  told  of  another  unknown  state.  This  fifth  unknown  cannot  be 
resolved  until  such  time  as  a  fifth  satellite  comes  into  view,  making  the  system 
solvable.  The  cycle  slip,  because  it  is  introducing  an  uncertainty  in  the  integer 


ambiguity  value,  causes  an  increase  in  the  overall  uncertainty  of  the  position.  The 
effect  of  the  slip  at  t=70  seconds  can  clearly  be  seen  in  the  single-frequency  case. 
The  uncertainty  curve  for  the  ideal  case,  that  without  tropospheric  and  ionospheric 
delay  biases,  is  not  effected  by  the  slip. 


FIGURE  8 

geometry  of  the  satellites 
relative  to  the  station 

Time  =  I  second 


FIGURE  9 

GEOMETRY  OF  THE  SATELLITES 
RELATIVE  TO  THE  STATION 
Time  =  60  seconds 


FIGURE  11 

GEOMETRY  OF  THE  SATELLITES 
RELATIVE  TO  THE  STATION 


Time  =  100  seconds 
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FIGURE  13 

SATELLITE  AZIMUTH  VERSUS  ELEVATION  OVER  TIME 
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FIGURE  14 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES  * 

Stationary  Satellite  and  Station 
Ideal  Case 
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FIGURE  15 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES 
Moving  Satellite  and  Station 
Ideal  Case 
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FIGURE  16 

UNCERTAINTY  USING  PSEUDO  RANGES 
w/  and  w/o  Tropospheric  and  Ionospheric  Errors 
Single-  and  Dual-Frequency  Models 
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FIGURE  17 

UNCERTAINTY  USING  PSEUDO  RANGES 
w/  and  w/o  Tropospheric  Errors 
Dual-Frequency  Model 
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FIGURE  18 

UNCERTAINTY  USING  PSEUDO  RANGES 
Ideal  Case  vs.  Real  World  Case 
Single-  and  Dual-Frequency  Models 
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FIGURE  19 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES 
w/  and  w/o  Tropospheric  and  Ionospheric  Errors 
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FIGURE  20 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES 
w/  and  w/o  Tropospheric  Errors 
Dual-Frequency  Model 
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FIGURE  22 

COMPARISON  OF  UNCERTAINTY 
USING  PSEUDO  AND  PHASE  RANGE  COMBINATIONS 


w/  and  w/o  Tropospheric  and  Ionospheric  Errors 
Ideal  Case  vs.  Real-World  Case 
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FIGURE  23 

UNCERTAINTY  USING  PSEUDO  RANGES 
Msmt  Model  w/  and  w/o  Elevation-Dependency 
Based  on  Troposphere 
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FIGURE  24 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES 
Msmt  Model  w/  and  w/o  Elevation-Dependency 
Based  on  Troposphere 
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FIGURE  25 

UNCERTAINTY  USING  PSEUDO  RANGES 
Msmt  Model  w/  and  w/o  Elevation-Dependency 
Based  on  Ionosphere 
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FIGURE  26 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES 
Msmt  Model  w/  and  w/o  Elevation-Dependency 
Based  on  Ionosphere 
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FIGURE  27 

UNCERTAINTY  USING  PSEUDO  RANGES 
Msmt  Model  w/  and  w/o  Radial  Orbit  Errors 
Coefficient  Error  =  10  meters 
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FIGURE  28 

UNCERTAINTY  USING  PSEUDO  RANGES 
Msmt  Model  w/  and  w/o  Radial  Orbit  Errors 


Coefficient  Error  =  1  meter 
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FIGURE  29 

UNCERTAINTY  USING  PSEUDO  RANGES 
w/  Elevation-Dependency  Based  on  Ionosphere 
w/  Tropospheric  and  Ionospheric  Errors 
Msmt  Model  w/  and  w/o  Radial  Orbit  Errors 
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FIGURE  30 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES 
Msmt  Model  w/  and  w/o  Radial  Orbit  Errors 
Coefficient  Error  =  10  meters 
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FIGURE  31 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES 
Msmt  Model  w/  and  w/o  Radial  Orbit  Errors 


Coefficient  Error  =  1  meter 
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FIGURE  32 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES 
w/  Elevation-Dependency  Based  on  Ionosphere 
w/  Tropospheric  and  Ionospheric  Errors 
Msmt  Model  w/  and  w/o  Radial  Orbit  Errors 
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FIGURE  33 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES  , 

Contribution  of  a  Single  Satellite's  Radial  Orbit  Errors 
w/  and  w/o  Tropospheric  and  Ionospheric  Errors 

► 


» 


) 


POSITION 

UNCERTAINTY 

(METERS) 


.0  10.0  40.0  40.0  00.0  100 


TIME  (SECONDS) 

FIGURE  34 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES 
Ideal  Case  vs.  Real-World  Noises 
Msmt  Noise(range)  =  Msmt  Noise(phase) 
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FIGURE  35 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES 
Msmt  Model  w/  and  w/o  Elevation-Dependency 
Based  on  Ionosphere 
Msmt  Noise(range)  =  Msmt  Noise(phase) 
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FIGURE  36 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES  1 

Msmt  Model  w/  and  w/o  Radial  Orbit  Errors 
Coefficient  Error  =10  meters 

Msmt  Noise(range)  =  Msmt  Noise(phase)  9 
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FIGURE  37 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES 
Msmt  Model  w/  and  w/o  Radial  Orbit  Errors 
Coefficient  Error  =  1  meter 
Msmt  Noise(range)  =  Msmt  Noise(phase) 
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FIGURE  38 

UNCERTAINTY  USING  PSEUDO  AND  PHASE  RANGES 
w/  Elevation-Dependency  Based  on  Ionosphere 
w/  Tropospheric  and  Ionospheric  Errors 
Msmt  Model  w/  and  w/o  Radial  Orbit  Errors 
Msmt  Noise(range)  =  Msmt  Noise(phase) 


CHAPTER  V 
CONCLUSIONS 


5.1  Overview 

This  thesis  examined  the  uses  of  different  types  of  GPS  data  together  with 
different  types  of  biases.  The  primary  theme  was  to  investigate  the  incorporation  of 
phase  range  measurements  into  pseudorange  measurements  with  the  use  of  the 
MSOFE  software  package.  Iterations  of  the  program  were  made  using 
combinations  of  a  variety  of  bias  sources  (integer  ambiguity,  tropospheric  delay, 
ionospheric  delay,  orbital  biases,  and  cycle  slips)  and  different  miscellaneous  models 
(equal  measurement  weightings,  elevation-dependent  measurement  noise  models). 
This  chapter  outlines  the  findings  of  this  research.  Areas  for  future  research  will 
also  be  addressed. 

In  this  research,  the  simulation  scenario  for  data  generation  utilized  a  specific 
broadcast  ephemeris,  station  location,  and  satellite  track.  Exact  results  will  be 
impossible  to  duplicate  under  different  satellite/station  geometries,  elevations,  and 
ephemeris  information.  The  findings  in  this  research  indicate  the  effects  of  the 
different  combinations  of  influences  and  may  be  used  to  surmise  the  effects  on 
scenarios  other  than  the  one  used  here. 


117 


5.2  Summary  jf  Results 


By  examining  the  results  of  the  research,  one  important  conclusion  is  clear: 
the  incorporation  of  phase  range  data  (known  to  a  millimeter)  along  with  the  meter- 
level  pseudorange  measurement  data  vastly  improved  the  overall  knowledge  of  the 
recovered  position  of  the  receiver.  With  the  addition  of  the  phase  range 
measurement,  the  uncertainty  of  the  position  decreased  (per  epoch)  on  the  order  of 
where  n  is  the  number  of  measurements 


The  addition  of  unmodeled  biases  in  the  "filter"  model  affects  one's  ability  to 
recover  the  station's  position.  In  all  cases,  the  bias-laden  models  generated  a 
position  uncertainty  greater  than  the  ideal  case  with  no  biases  present  Biases  in  the 
troposphere  did  not  increase  the  uncertainty  level  over  the  ideal  case  by  as  great  a 
magnitude  as  the  ionospheric  bias  did,  because  the  tropospheric  delay  was  known  to 
98%  and  the  ionospheric  delay  only  to  50%. 

The  elevation-dependent  model,  based  on  the  troposphere,  showed  little 
change  in  the  uncertainty  when  only  pseudoranges  were  used.  The  (troposphere) 
elevation-dependent  model  did,  however,  indicate  an  improvement  over  the 
elevation-independent  mode!  for  the  pseudo/phase  range  combinations.  The  large 
rise  in  position  uncertainty  in  the  latter  model  due  to  the  redundant  satellite  was 
absent  in  the  former  model,  due  to  the  adequate  weighting  of  the  low  (in  elevation) 
satellites.  For  the  ionosphere  elevation-dependent  model,  the  pseudorange-only 
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measurement  model  did  not  see  improvement  due  to  the  redundant  satellite.  The 
uncertainty  actually  rose  due  to  the  trade-off  between  the  weighting  of  the  low 
satellites  and  the  unfavorable  geometry  caused  by  essentially  omitting  the  fifth 
measurement.  The  pseudo/phase  range  combination,  again,  showed  an 
improvement  over  both  the  pseudorange-only  model  and  the  pseudo/phase  range 
model  without  the  elevation-dependency.  However,  the  degrading  geometry  due  to 
omitting  the  redundant  satellite  could  cause  a  problem  in  the  position  recovery. 

The  presence  of  the  orbital  biases  served  to  increase  the  uncertainty  in  both 
measurement  models.  The  greater  the  a  priori  uncertainty  in  the  orbital  bias 
coefficients,  the  greater  the  uncertainty  in  the  receiver's  position.  As  the  a  priori 
uncertainty  of  the  bias  increased  to  + 1  m,  the  overall  uncertainty  curve  from  the 
filter  approached  that  of  the  single-frequency  model  which  was  void  of  orbital 
biases. 

When  the  measurement  noises  of  the  pseudoranges  and  phase  ranges  were 
equal,  the  station's  position  uncertainty  for  the  five-satellite  epochs  was  improved 
over  the  unequal-noise  case  (with  presence  of  no  other  biases).  However,  as  the 
geometry  degraded  due  to  the  movement  of  the  satellites,  the  equal-noise  case  could 
not  recover  the  position  to  an  uncertainty  as  good  as  the  unequal-noise  case.  The 
latter  case,  using  the  phase  ranges  good  to  +1  mm,  produced  better  results  in 
epochs  with  poor  satellite/station  geometries. 
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When  orbital  errors  were  incorporated  into  the  equal-measurement-noise 
model,  the  uncertainty  curves  resembled  those  of  pseudorange-only  models.  The 
worst  uncertainty  of  all  was  discovered  in  the  single-frequency,  equal-measurement- 
noise,  elevation-dependent  (based  on  ionosphere)  case.  The  uncertainty  approached 
1 26  m  at  its  peak  and  showed  no  signs  of  improvement  as  the  redundant  satellite 
came  into  view. 

Overall,  the  incorporation  of  mm-level  phase  range  measurements  into  the  m- 
level  pseudorange  measurements  served  to  reduce  the  uncertainty  in  solving  for  the 
unknown  variables  in  the  state  vector.  As  more  corruption  was  introduced  into  the 
measurement  model,  the  ability  to  recover  accurately  the  receiver's  position  was 
degraded.  Often,  the  elevation-dependent  measurement  noise  model  further  reduced 
the  uncertainty  of  the  station  position,  providing  that  the  satellite/station  geometry 
created  from  omitting  a  measurement  was  sufficient  enough  for  adequate  recovery. 
Certainly  a  trade-off  exists  in  reducing  the  weighting  of  a  measurement  if  the 
geometry  of  the  remaining  satellites  cannot  support  the  rest  of  the  system.  This 
research  set  the  elevation  cut-off  at  the  horizon  (0  degrees).  Better  results  could  be 
achieved  if  the  cut-off  were  raised  to  a  greater  level,  e  g.,  15  degrees.  The 
ionospheric  and  tropospheric  delays  at  0  degrees  were  often  prohibitive  in  one's 
ability  to  recover  positions  to  an  adequate  level. 
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5.3  Areas  for  Future  Research 

This  thesis  concentrated  on  the  single-frequency  GPS  receiver  to  model  the 
measurements  The  dual-frequency  receiver  may  serve  to  provide  better  results.  In 
the  dual-frequency  receiver,  biases  such  as  the  ionosphere  (with  a  +50% 
uncertainty)  could  be  eliminated  from  the  problem  altogether  by  combining  in  a 
particular  way  the  two  measurements  at  the  respective  frequencies.  This 
combination,  however,  has  the  drawback  of  amplifying  measurement  noise.  The 
resulting  noise  value  is  scaled  up  by  approximately  3  times  the  original  L,  or  L, 
noise  value*.  The  trade-off  may  well  be  worth  the  added  noise,  because  the 
ionospheric  delay  is  very  difficult  to  model,  and  all  attempts  to  model  the 
ionospheric  effects  will  probably  result  in  errors  very  much  greater  than  a  three-fold 
increase  in  measurement  noise. 

Further  research  in  this  area  also  could  include  the  incorporation  of  Doppler 
measurements  (time  rate  of  change  of  phase  (ranges))  into  the  filter  process 
augmenting  the  pseudoranges  and  phase  ranges.  As  seen  in  this  research,  the 
inclusion  of  better-known  information  serves  well  to  increase  the  accuracy  of  the 
recovered  information.  An  added  plus  to  including  Doppler  measurements  is  that 
there  are  no  integer  ambiguity  problems  to  contend  with.  A  trade-off  must  occur. 


'The  noise  on  the  L,  carrier  is  amplified  by  2.5  times,  and  the  noise  on  the  L;  is  amplified  by  1.5 
times.  The  overall  effect  on  the  standard  deviation,  a,  is 

TOm,/l.52+2.52=  2.91a 


a  =a 
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however.  Using  Doppler  measurements,  one  must  attempt  to  solve  for  the  rate  of 
change  of  clock  drift  (i.e.,  another  unknown  must  be  added  to  the  state  vector  to 
account  for  the  clock).  One  might  "solve"  this  situation  by  increasing  the  precision 
of  the  clock  to,  perhaps,  that  of  a  rubidium  clock,  known  to  1  part  in  1012.  This 
way,  the  unknown  frequency  drift  of  the  clocks  could  be  "controlled"  and  act  as  a 
stabilizing  influence  in  recovering  the  position  of  the  receiver. 

One  final  issue  must  be  addressed— the  issue  of  an  updated  Kalman  gain  matrix 
for  the  augmented,  "real-world"  situation.  Recall  from  Chapter  2  that,  with  and 
without  the  presence  of  real-world  errors,  the  covariance  update  equations  take  on 
different  forms  (see  Equations  2-45  and  2-39,  respectively).  The  program  MSOFE, 
however,  did  not  update  the  Kalman  gain  matrix  (Equation  2-40)  to  account  for  the 
additional  vector  of  errors.  In  all  cases,  the  Kalman  gain  used  by  MSOFE  in  the 
update  process  was  the  usual  filter  gain,  shown  in  Equation  4-11. 

K  =  P'(H~)t[H'P"(H~)t  +  R]~'  (5-1) 

However,  this  value  did  not  take  into  account  any  of  the  additional  real-world 
error  states.  Solving  for  the  Kalman  gain  matrix  based  on  augmented  predicted 
residual  variances,  one  gets 

K augmented  =  P' (H‘  )T[H~p- (H' )T  +  R  +  MPyMT ]  ' 1  (5-2) 


which  incorporates  the  covariance  of  the  error  states  into  the  update  process. 
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Further  study  should  examine  the  accuracy  of  the  results  using  MSOFE 
together  with  this  augmented  gain  matrix.  In  MSOFE,  the  Kalman  gain  matrix  in 
Equation  4-11  was  employed.  The  R  matrix  was  a  diagonal  matrix  due  to  the  lack 

of  correlation  of  noise  between  states.  By  adding  the  vector  of  "real-world"  errors, 
Y ,  to  the  system,  one  gets  an  augmented  measurement  noise,  seen  in  Equation  4-13. 

Ra„d=R  +  MPyMT  (5-3) 

This  matrix  is  clearly  no  longer  diagonal.  Here,  the  noise  becomes  correlated 
between  states— a  condition  ignored  by  MSOFE,  which  assumed  the  same  systematic 
effects  for  all  states.  If  the  same  systematic  effects  did  occur  epoch  after  epoch,  one 
could  investigate  including  these  y  states  in  the  "filter"  model  to  attempt  to  balance 

the  "system"  and  "filter"  models  properly  in  order  to  gain  better  information  about 
them.  Although  MSOFE's  methods  were  satisfactory  for  investigating  the  topics  in 
this  thesis,  they  were  by  no  means  optimum. 

The  issue  of  the  Pr  matrix  should  also  be  addressed.  In  Chapter  2,  we 

assumed  that  this  covariance  matrix  remained  bounded.  However,  one  should  not 
ignore  the  situation  where  some  states  in  Pr  become  unbounded.  As  mentioned 
above,  one  may  desire  to  treat  some  states  in  the  y  vector  as  unknowns  and  attempt 

to  solve  for  them2.  This  situation  may  set  up  a  condition  whereby  the  covariance 
values  increase  rapidly  (i.e.,  PY-»oo).  In  this  case,  the  augmented  R  matrix  from 

2  Techniques  concerning  the  elimination  of  parameters  may  be  found  in  the  work  bv  Prijatna, 

1992. 
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Equation  5-3  will  also  become  unbounded  (for  those  states),  and  the  Kalman  gain 
matrix  will  approach  zero  (for  those  states).  One  cannot  easily  predict  the  exact 
effect  on  R  and  K  that  an  unbounded  covariance  matrix  will  have;  one  may  or  may 
not  be  able  to  solve  for  the  desired  states  in  the  y  vector.  Further,  if  the  gain  matrix 

approaches  zero,  then  the  states  will  not  be  updated  from  the  current  states  (i.e.,  one 
would  not  see  a  reduction  in  P  from  P"  toP+).  Further  research  in  this  area  could 
investigate  the  effect  of  this  situation  on  the  outcome  of  the  uncertainty  of  the 
unknowns  and  determine  if  this  situation  is  detrimental  to  the  adjustment  procedure 
and  the  recovery  of  accurate  information. 
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